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ABSTRACT 


This  study  is  directed  at  the  problem  of  identifying  the 
processes  involved  in  causing  deflagrating  propellants  to  deton- 
ate, and  use  the  results  to  devise  means  for  minimizing  rocket 
motor  detonations  consistent  with  optimal  motor  performance  and 
reliability.  Of  primary  concern  to  this  study  are  rapid  pressure 
rises  within  cracks,  flaws  or  debonds  within  propellant  grains 
that  could  Initiate  detonation. 

This  report  covers  the  formulation  of  a computerized  model 
with  which  to  study  the  consequence  of  sudden  exposure  of  cracks 
to  gas  at  higher  pressures  and  temperatures  than  those  initially 
existing  within  the  crack.  Simple  one-dimensional  formulations 
were  developed  to  predict  the  transient  gas  flows  and  propellant 
regression  rates  within  cracks  of  variable  shapes  with  the  down- 
stream end  of  the  crack  open  or  closed.  Provisions  were  also 
made  for  predicting  the  consequence  of  reflected  stress  waves  upon 
the  crack  thicknesses. 

The  most  important  finding  of  this  study  is  that  stress  waves 
reflected  from  stiff  motor  cases  can  cause  cracks  to  contract  at 
sufficient  rates  to  produce  appreciable  pressure  rises.  Unfor- 
tunately the  maximum  pressures  and  pressure  rises  could  not  be 
quantified  due  to  computational  difficulties  in  the  gas  dynamic 
calculations  created  by  highly  variable  crack  thicknesses  produced 
during  the  final  stages  of  crack  contraction.  Questions  regarding 
possible  impact  of  burning  crack  surfaces  remain  to  be  resolved. 
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1.  INTRODUCTION 


During  recent  years  high-energy  propellants  have  been  devel- 
oped containing  HMX  in  order  to  improve  the  performance  of  rocket 
motors.  During  the  course  of  this  development  two  explosions  have 
occurred.  The  first  was  the  Hercules  plant  in  Bacchus,  Utah  in 
May  1974,  the  second  was  at  the  Naval  Weapons  Center  in  China  Lake, 
California.  Experimental  evidence  (1,2,3)  indicates  that  a first 
necessary  condition  for  detonation  is  a rapid  pressure  rise.  For 
this  reason,  this  report  is  concerned  with  identifying  the  condi- 
tions by  which  burning  within  cracks,  flaws  or  debonds  will  result 
in  pronounced  pressure  rises. 

Deflagration  to  detonation  (DDT)  can  occur  as  a result  of  a 
number  of  stimuli  such  as  heat  and  impact  (1,4-6).  Such  stimuli 
can  cause  a burning  reaction  that  rapidly  accelerates  under  suit- 
able conditions  causing  shock/stress  wave  formation  and  subsequent 
detonation,  The  existing  detonation  theories  describe  the  initia- 
tion of  propellants  and  explosives  by  strong  shock  waves  reasonably 
well.  However,  detonation  can  occur  as  a result  of  low-amplitude 
shocks  that  are  well  below  those  predicted  by  hydrodynamic  theories 
of  detonation.  Equally  important  as  high-order  detonations  is 
the  class  of  chemical  decompositions  referred  to  as  subdetonation 
reactions.  While  such  reactions  are  not  as  destructive  as  high- 
order  detonations,  they  can  destroy  rocket  motors. 

In  view  of  the  complex  nature  of  this  problem,  a combined 
theoretical/experimental  investigation  is  called  for.  Analytical 
studies  reveal  parameters  whose  role  or  significance  is  frequently 
lost  or  not  measurable  in  experiments;  while  experiments  serve 
to  validate  the  accuracy  of  analytical  models  and  to  resolve  any 
untractable  features  of  the  problem. 

Recognizing  the  importance  of  putting  the  problem  into  per- 
spective, we  have  focused  our  attention  upon  the  development  of 
simple  analytical  models  with  which  to  gain  a better  appreciation 
of  the  importance  of  various  parameters  in  causing  the  initiation 
of  high-energy  propellants.  Through  such  studies  a clearer  pic- 
ture will  be  gained  by  which  to  identify  the  kinds  of  analysis 
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and  experiments  needed  to  compliment  and  guide  model  development, 
and  to  speed  the  resolution  of  this  complex  problem. 

Important  aspects  of  the  problem  requiring  consideration  are 
transient  heating  and  burning  of  the  propellant,  gas  dynamics, 
crack  propagation,  and  positive  and  negative  deformations  of  cracks 
that  retard  or  speed  burning.  For  this  purpose  computer  routines 
have  been  developed  using  a number  of  simplifying  assumptions  and 
datum  estimates  --  none  of  which  is  believed  will  compromise  iden- 
tification of  the  essential  features  of  the  problem.  Particular 
attention  was  given  to  structuring  the  computerized  model  so  it 
can  be  readily  upgraded  --  serving  as  it  does  as  a research  tool. 

The  most  important  findings  of  this  study  are  twofold.  The 
first  is  the  consequence  of  sensible  heat  deposited  in  the  propel- 
lant during  the  relatively  long  times  prior  to  appreciable  pressure 
rises.  Such  heat  will  speed  burning  over  and  above  quasi-steady 
state  values  if  and  when  thermal  fluxes  increase  due  to  rising  gas 
pressures,  velocities  and/or  temperatures. 

Of  equal  or  greater  importance  is  closure  of  cracks  by  stress 
waves  reflected  from  motor  cases  or  generated  by  other  cracks . 

Under  suitable  conditions,  i.e.  stiff  case  and  sudden  exposure 
of  cracks  to  high  pressure  gas,  cracks  will  contract  causing 
pronounced  pressure  rises. 

One  other  possibility,  which  has  not  been  explored  in  this 
study,  is  the  fragmentation  of  the  propellant  by  multiple  crack- 
ing resulting  in  burning  fragments  being  propelled  against  one 
another.  In  this  regard  high  explosives  can  be  fragmented  during 
high-pressure  burning.  Evidence  of  this  is  shown  in  Fig.  1 from 
Ref.  8. 


Fragments  shown  in  Fig.  1 were  recovered  from  a closed  bomb 
containing  a burning  cylinder  of  Composition  B following  separa- 
tion of  the  closed  bomb.  Failure  was  caused  by  abnormal  in- 
creases of  pressure,  which  in  other  similar  teats  lead  to  detona- 
tion (1).  In  this  test,  burning  was  extinguished  by  the  sudden 
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drop  in  pressure.  Ir  was  concluded  that  all  surfaces  of  all 
Composition  B fragments  were  in  the  process  of  burning  as  evi- 
denced by  a "frozen"  foam  layer  over  all  surfaces.  In  other 
similar  tests  (1)  it  was  observed  that  ionization  probes  located 
at  various  depths  within  the  explosive  fired  in  a random  fashion 
with  respect  to  depth.  Random  firing  occurred  during  the  period 
of  abnormal  pressure  rises  starting  at  about  1000  psi,  and  fre- 
quently was  followed  by  detonation.  One  possible  explanation 
of  the  random  firing  is  multiple  cracking  such  as  indicated  by 
Fig.  1. 

In  the  remainder  of  this  report  we  shall  discuss  analyses 
used  in  developing  the  model  and  results  produced  thereby.  At 
the  conclusion  we  shall  indicate  which  of  the  many  areas  of  un- 
certainties are  in  greatest  need  of  study. 
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2.  MODEL  DEVELOPMENT 


Before  commencing  model  development,  the  literature  was  re- 
viewed for  information  and  data  to  assist  in  describing  various 
aspects  of  the  DDT  process,  and  to  guide  the  conduct  of  the  pro- 
gram, In  this  regard,  considerable  information  was  found  describ- 
ing specific  aspects  of  the  problem.  This  includes 

• Experimental  ignition  of  propellants  and 
explosives  by  shocks,  radiation  , hot 
wires  and  hot  gases  (9-11) 

• Experiments  describing  the  spread  of 
burning  or  flame  over  surfaces  or  into 
pores  or  cracks  (12-14) 

• Combustion,  ignition  and  detonation 
theories  and/or  experiments  (15-30) 

• Studies  of  crack  initiation  and  propaga- 
tion (1,31,32) 

Unfortunately  no  comprehensive  treatment  of  the  problem  of  crack 
. deformation  and  burning  was  found  with  which  to  put  the  problem 
into  perspective.  This  is  not  at  all  surprising  in  view  of  the 
many  complex  interdependent  phenomena  involved. 

Based  upon  this  finding  a dynamic  model  was  developed  using 
a number  of  simplifications  to  gain  a better  appreciation  of  the 
problem.  In  this  model  provisions  were  made  for 

• Single  crack  within  monopropellant  of  variable 
dimensions,  shape  and  end  conditions 

• Sudden  flow  of  hot  high-pressure  gases  (with 
friction)  into  crack 

• Crack  or  flaw  propagation  (at  present  intro- 
duced artificially) 

• Crack  deformations  --  positive  and  negative 
--  induced  by  pressure  environment  and  case 
response  to  stress  waves 


• Convective  (turbulent)  and  internal  heating  of 
propellant 

« Propellant  regression  rates  due  to  one  of  two 
idealized  situations  namely 

- gas  evolution  assuming  molten  propellant 
remains  in  place 

- immediate  discharge  of  molten  propellant 
into  gas  stream 

At  present,  gas  flow  is  treated  as  one-dimensiona  /ith  friction, 
assuming  that  the  gas  behaves  as  a perfect  gas  and  that  reactions 
occur  instantaneously.  Propellant  temperatures  are  calculated 
on  a transient  basis  and  used  to  predict  regression  rates  at  vari- 
ous crack  locations  for  either  one  of  the  two  idealized  caseB 
cited  above.  The  consequences  of  transient  pressures  upon  crack 
deformation  is  modeled  assuming  the  propellant  behaves  as  a simple 
linear  isotropic  visco-elastic  material.  Yield  or  failure  of 
surrounding  propellant  is  not  considered,  and  the  dependence  of 
property  data  upon  temperature  and  pressure  is  neglected.  More- 
over, no  provisions  have  been  made  for  heat  feedback  from  the 
burning  propellant  other  than  through  the  effect  of  burning  upon 
the  bulk  presaures,  temperatures  and  velocities  of  the  gases. 

Figure  2 illustrates  provisions  made  regarding  the  crack 
and  flow  conditions.  The  crack  was  subdivided  into  uniformly 
long  elements.  -The  width  of  the  crack  is  considered'  sufficient 
to  accomodate  a one -dimensional  treatment  of  the  gas  flow.  Means 
were  provided  for  variable  crack  thicknesses  from  element  to 
element.  In  addition  provisions  were  provided  for  altering  the 

e Initial  pressure  and  temperature  within  the  crack 

e End  conditions  of  crack  --  open  or  closed 

e Gas  pressure  and  temperature  in  cavity  at  the 
upstream  end  of  the  crack 

e Length  of  crack  as  a function  of  time  (crack 
propagation) 

In  tha  remainder  of  this  section  we  shall  discuss  how  each 
of  the  principal  phenomena  has  been  modeled  and  pertinent  results 
obtained  thereby.  Topics  covsrad  are  gas  dynamics,  mechanical 
response  of  cracks  and  transient  regression  velocities  of  propellant. 
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Figure  2 DESCRIPTION  OF  CRACK  CONDITIONS  CONSIDERED 


2 . 1 Gas  Dynamics 

2.1.1  Gas  Dynamic  Model 


One  essential  aspect  of  the  problem  is  the  transfer  of  gas- 
eous products  and  heat  along  the  length  of  the  crack.  To  achieve 
high  convective  heating  and  strong  stress-waves,  we  have  chosen 
to  consider  the  problem  of  suddenly  exposing  a crack  to  a high 
pressure  gas  cavity. 

Even  though  the  gas  dynamic  model  includes  a variety  of  im- 
portant flow  features,  it  represents  a preliminary  version  subject 
to  modification.  These  modifications  will  be  identified  and  in- 
corporated in  the  course  of  treating  various  aspects  of  the  problem. 

The  gas  dynamic  model  is  based  upon  an  Eulerian  representa- 
tion of  the  flow  region  and  is  limited  to  a one-dimensional  treat- 
ment. The  selection  of  the  Eulerian  representation  permits  a 
simple  identification  with  locations  along  the  crack  and  allows 
for  mixing  of  freshly  burnt  reaction  gas.  It  is  also  convenient 
for  large  movements  experienced  by  some  of  the  gas. 

The  one-dimensional  treatment  is  appropriate  because  the 
width  and  length  of  cracks  are  usually  large  compared  to  their 
thickness.  Thus  the  most  significant  flow  gradients  and  move- 
ments will  occur  along  the  length  of  the  crack.  Since  flow 
gradients  will  exist  across  the  thickness  of  the  crack,  the 
values  of  the  variables  used  in  the  one-dimensional  treatment 
represent  average  cross-sectional  values. 

The  model  treats  a crack  of  Borne  specified  thickness  varia- 
tion at  zero  time.  The  gas  in  the  crack  at  zero  time  is  con- 
sidered to  have  the  same  composition  as  the  reaction  products, 
and  to  exist  in  a rest  state  (no  flow)  at  some  specified  pres- 
sure and  temperature.  Tha  crack  is  then  suddenly  connected  to  a 
cavity  filled  with  gaseous  reaction  products  at  high  temperature 
and  pressure.  The  cavity  state  is  treated  as  a constant  with 
respect  to  time  due  to  the  relatively  short  duration  of  the  events 
occurring  within  the  crack. 
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The  downstream  end  of  the  crack  can  either  be  closed  or  con- 
nected to  another  cavity  held  at  a specified  low  pressure.  The 
low-pressure  cavity  could  be  a large  debonded  region.  For  open- 
ended  cracks,  the  crack  length  is  held  constant.  If  the  closed 
crack  condition  is  selected,  the  crack  may  be  allowed  to  grow 
in  length  in  some  arbitrary  fashion.  At  present  the  growth  fea- 
ture is  not  based  upon  any  physically  plausible  crack  propagation 
mechanism.  The  incorporation  of  such  a criterion  together  with 
crack  branching  represents  one  of  the  ultimate  goals  for  improving 
the  model. 

2.1.2  Gas  Dynamic  Results 

A computer  program  was  developed  to  implement  the  gas  dynamic 
model.  This  computer  program  exists 

(1)  in  the  form  of  a main  program  with  which  to 
examine  some  aspects  of  the  gas  dynamic  phen- 
omena and  to  Berve  as  a vehicle  for  further 
development,  and 

(2)  in  the  form  of  a subroutine  of  the  overall 
Dynamic  Model  of  Crack  Burning. 

These  codes  were  programmed  in  FORTRAN  IV  user-oriented  source 
language  and  have  been  successfully  executed  on  a UNIVAC  1108, 
a large  scientific  computer, 

This  section  presents  the  computed  results  of  a reference 
problem  to  acquaint  the  reader  with  the  adequacy  of  the  model  and 
with  a typical  response  of  the  gas  dynamic  system  in  the  absence 
of  a reacting  propellant.  The  reference  problem  was  selected 
for  the  validation  of  the  code  and  the  subsequent  evaluation  of 
various  effects.  The  problem  configuration  is  that  of  a constant 
thickness  crack,  which  is  10  cm  long  with  the  downstream  end 
closed.  This  crack  is  connected  at  the  upstream  end  to  a cavity 
filled  with  air  at  a pressure  of  60  bars  and  a temperature  of 
2000 °C.  Air  is  also  in  the  crack  or  flow  channel,  but  is  at  a 
pressure  of  1 Bar  and  a temperature  of  20*0.  The  thickness  of 
the  crack  will  not  be  a factor  until  certain  physical  effects 
such  a friction  are  introduced.  The  initial  validation  examines 


the  inertial  effects  of  the  flow.  The  influence  of  the  various 
other  effects,  such  as  friction,  area  variation,  etc.  have  been 
examined,  but  are  not  diocussed  herein. 

At  time  equal  to  zero  a shock  wave  will  be  generated  at  the 
cavity  end  of  the  flow  channel  and  propagate  downstream  compres- 
sing the  gas  which  was  originally  in  the  crack  to  a pressure  of 
20.4  Bars.  The  cavity  gas  will  expand  to  a sonic  state  at  the 
inlet  region  of  the  crack  and  then  undergo  a further  (nonsteady) 
expansion  within  the  crack  until  the  pressure  drops  to  20.4  Bars. 
Thus  a centered  rarefaction  wave  system  will  also  be  generated 
at  the  inlet  and  propagate  downstream.  The  corresponding  wave 
diagram  is  presented  in  Figure  3.  An  interface  moving  along 
with  the  fluid  separates  the  two  masses  of  gas.  Region  A is 
the  undisturbed  region;  Region  B is  the  compressed  gas  state  re- 
gion; Region  C is  the  fully  expanded  cavity  gas  region;  and  Region 
D is  the  nonsteady  expansion  region.  Regions  A,  B,  and  C are 
uniform  state  -regions. 

The  shock  wave  is  travelling  at  a velocity  of  139.8  cm/ms 
and  accelerates  the  rest  gas  to  a velocity  of  109.9  cm/ms.  The 
fully  expanded  cavity  gas  is  also  accelerated  to  this  velocity. 
These  two  uniform  regions  are  at  the  same  pressure  and  flow  velo- 
city conditions;  however,  their  entropy  is  different.  The  tem- 
perature of  the  compressed  gas  is  1007°C  corresponding  to  a sound 
velocity  of  69.6  cm/ms.  The  cavity  gas  sound  velocity  is  92.7 
cm/ms  in  the  cavity,  84.6  cm/ms  at  the  sonic  inlet  state  and  79.6 
cm/ms  in  the  fully  expanded  state  (Region  C) . The  gas  pressure 
at  the  sonic  inlet  state  is  31.7  Bars.  This  initial  flow  field 
is  self-similar  in  nature  for  a short  period  of  time  in  that  all 
the  variables  are  constant  along  rays  passing  through  the  origin. 
This  period  of  time  terminates  when  the  initial  Bhock  wave  re- 
flects from  the  closed  end  of  the  flow  channel  and  creates  a 
uniform  high  pressure  rest  region,  Region  E.  The  reflected  shock 
wave  propagates  back  upstream  and  interacts  with  the  Interface. 


This  shock-interface  interaction  generates  additional  disturbances 
and  a more  complex  flow  environment  is  established. 

The  initial  self-similar  flow  field  can  be  used  to  evaluate 
the  fidelity  of  the  computer  flows  during  the  entry  phase.  The 
lengthwise  or  distance  resolution  will  depend  upon  the  number  of 
cells,  MT,  used  to  define  the  flow  channel.  In  general  20  cells 
have  been  used  to  define  the  10  cm  long  crack,  although  in  some 
preliminary  calculations  only  10  cells  have  been  used.  Thus,  for 
the  former  case  the  cells  are  0.5  cm  long.  For  a given  crack 
length  the  number  of  calculation  units  required  to  obtain  a solu- 
tion for  a given  period  of  time  will  be  proportional  to  MT^  since 
MT  effects  both  the  distance  and  time  increments. 

The  pressure  and  velocity  distribution  within  the  crack  are 
presented  in  Figure  4.  The  dotted  line  represents  the  theoretical 
solution  at  the  corresponding  time.  These  resultB  are  for  Time 
Step  20  and  corresponds  to  a time  of  0.044  ms.  Approximately  12 
of  the  20  cells  have  been  engulfed  by  the  shock  wave  at  this  time, 
hence  this  comparison  represents  a rather  severe  test  of  the 
fidelity  of  the  calculations.  Notice  that  the  shock  discontinuity 
(located  at  a distance  of  6.2  cm)  is  smeared  over  some  small  dis- 
tance. Similar  results  of  the  pressure  distribution  In  the  vicinity 
of  the  shock  front  are  shown  in  the  inset  for  the  case  where  40 
cells  have  been  used  to  define  the  crack.  These  latter  results 
are  better  at  the  shock  front  and  support  the  known  fact  that 
shock  waves  will  be  smeared  out  over  three  to  four  cells  when 
this  type  of  numerical  method  is  used.  There  is  also  some  over- 
shoot in  the  vicinity  of  the  shock  front.  The  degree  of  overshoot 
varies  with  the  flow  variables.  The  front  (downstream  edge)  of 
the  centered  rarefaction  wave  is  located  at  a distance  of  1.3  cm. 

The  computed  results  for  this  flow  feature  are  more  diffuse  than 
the  theoretical  solution,  however  this  should  be  expected  since 
only  several  ceils  are  used  to  define  this  wave  system  detail  at 
this  early  time. 
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The  pressure  and  particle  velocity  are  correspondingly  equal 
in  Regions  B and  C;  however,  the  other  variables  such  as  tempera- 
ture, density,  etc.  are  discontinuous  at  the  interface.  The  in- 
terface at  this  time  is  located  at  4.8  cm.  The  numerical  results 
reproduce  these  discontinuities  with  about  the  same  fidelity  as 
they  reproduced  the  shock  discontinuity. 

The  reflected  shock  wave  propagates  back  up  toward  the  high 
pressure  cavity  end  of  the  crack  and  interacts  with  both  the 
interface  and  the  centered  rarefaction  wave.  The  gross  effect 
of  this  reflection  process  is  to  decelerate  the  gas  to  essentially 
a rest  state  and  to  increase  the  pressure  to  approximately  105 
Bars.  This  strong  disturbance  will  Interact  with  the  upstream 
end  of  the  crack  and  create  a more,  diffuse  rarefaction  wave  system 
which  will  then  propagate  downstream  initiating  another  cycle  of 
the  transient  flow.  This  second  cycle,  however  represents  a 
pressure  drop  as  compared  to  the  pressure  rise  (shock  wave)  of 
the  first  cycle.  The  primary  disturbance  will  continue  to  rever- 
berate within  the  crack  producing  alternate  cycles  of  compression 
waves  and  rarefaction  waves  of  decreasing  magnitude  until  a ques- 
cent  or  rest  state  is  achieved  which  is  in  pressure  equilibrium 
with  the  high  pressure  cavity.  Temperature  and  density  gradients 
will  exist  within  the  crack  due  to  the  nonuniform  entropy  produc- 
tion, however  these  variations  will  be  rather  small.  The  time 
history  of  the  pressure,  particle  velocity  and  temperature  at  the 
two  ends  of  the  crack  are  presented  in  Figures  5,  6,  and  7 res- 
pectively. The  particle  velocity  at  the  closed  end  is  Identically 
zero  and  is  therefore  not  shown. 

2.2  Mechanical  Response  of  Cracks  (Lateral) 

2.2.1  Mechanical  Response  Model 

The  flow  of  high  pressure  gases  into  cracks  within  and  adja- 
cent ot  the  propellant  mass  and  any  subsequent  burning  of  these 
propellant  surfaces  will  create  a time  varying  pressure  environ- 
ment with  the  cracks.  The  surrounding  propellant  mass,  which 
was  originally  in  mechanical  equilibrium  with  a low  pressure  dis- 
tribution with  the  crack  will  respond  to  these  new  mechanical  loads. 
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Figure  6 PARTICLE  VELOCITY  HISTORY  - REFERENCE  SOLUTION 
CELL  2 


Figure  7 TEMPERATURE  HISTORIES  - REFERENCE  SOLUTION 


Two  mechanical  response  models  have  been  developed  and  are 
used  in  conjunction  with  the  gas  dynamic  mode.  One  response  model 
treats  the  early  phases  of  the  mechanical  response  when  stress 
waves  radiate  away  from  the  crack  and  the  crack  is  in  a growth 
phase.  The  second  model  treats , in  a simple  fashion,  the  inter- 
action of  these  radiating  waves  with  the  confining  shell  structure 
and  their  subsequent  interaction  with  the  crack.  Under  these  con- 
ditions the  crack  growth  may  be  arrested  and  lead  to  partial  or 
complete  crack  collapse. 

2.2.2  Mechanical  Response  Results 

The  mechanical  response  of  the  propellant  mass  to  the  pres- 
sure environment  within  the  crack  is  very  complex  due  to  the 
multiplicity  of  stress  wave  types  present  and  the  variety  and 
complexity  of  the  overall  geometry  of  the  system,  i.e. , the  crack, 
the  propellant  mass,  and  the  case  configuration  and  location. 

The  stress  wave  system  in  the  propellant  mass  during  the  initial 
transient  is  illustrated  in  Figure  8.  The  primary  motion  of 
the  propellant  will  be  in  a direction  normal  to  the  crack  and 
will  be  roughly  one  of  plane  strain,  For  this  reason  a simple 
local  plane  strain  representation  for  the  transient  stress  field 
was  adopted  as  a model  of  the  mechanical  response  of  the  propel- 
lant mass  in  the  absence  of  remote  boundary  effects  such  as  the 
case  response.  In  this  manner  the  crack  thickness  growth  feature 
could  be  coupled  directly  into  the  gas  dynamic  model  and  the 
coupled  response  determined,  ThiB  initial  response  model  is 
valid  until  the  stress  waves  which  radiate  away  from  the  crack 
Interact  with  the  remcw,e  boundaries  of  the  propellant  mass  and 
its  confining  system,  and  propagate  back  to  the  crack  locations. 
This  process  is  very  complex  and  varied. 

A simple  shell  response  model  was  selected  as  an  initial 
version  of  the  remote  boundary  effect,  The  model  treats  the 
shell  as  a simple  single  degrec-of-freedom  lumped  mass  which 
is  position  restreined  and  which  responds  to  the  mean  radiative 
stress  field.  The  response  generates  a reflected  streBs  wave 
syatem  which  arrives  beck  at  the  crack  after  some  specified 
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Traveling  Pressure  Pulse 


Figure  3 STRESS  WAVE  SYSTEM  IN  PR< 


delay  period  corresponding  to  the  mean  distance  to  the  remote 
boundary.  This  model  is  illustrated  in  Figure  9.  Two  typical 
crack  thickness  histories  are  shown  in  Figure  10  for  a cell 
near  the  midpoint  of  the  crack  for  (1)  a cased  propellant  mass 
and  (2)  one  where  the  remote  boundary  is  very  far  away  (i.e., 
a very  large  propellant  mass)  from  the  crack. 

2 . 3 Dynamic  Burn  Model 

2.3.1  Problem  Definition 

Before  discussing  the  dynamic  model  it  is  desirable  to  ex- 
amine two  features  of  the  problem  having  to  do  with  crack  sur- 
faces. The  first  concerns  the  amount  of  exposed  surface  area, 
and  the  second  the  limitations  of  a one-dimensional  treatment. 

Exposed  Surface  Area 

Crack  surfaces  are,  of  course,  highly  irregular  and  will 
present  greater  areas  than  that  of  a planar  surface.  To  appre-  , 
ciate  the  extent  by  which  crack  areas  exceed  planar  areas,  we 
shall  consider  the  two  crack  surfaces  to  contain  a number  of  spher 
ical  particles  of  various  radii  r.  The  particles  shall  have  half 
or  less  of  their  area  exposed.  In  that  the  exposed  area  is  propor- 
tional to  the  height  x of  the  particle  above  base  plane  associated 
with  the  crack,  the  exposed  particle  area  is  given  by  2irrx 
for  either  crack  surface.  On  the  other  hand  the  cross-sectional 
area  of  particles  of  radius  r in  the  base  plane  is  given  by 
tr(2rx-x  ).  Since  all  x values  from  0 to  r are  equally  likely, 
the  average  areas  associated  with  each  particle  can  be  found  by 
integrating  each  of  the  above  expressions  with  respect  to  x and 
dividing  by  r.  The  ratio  of  the  average  areas  of  the  exposed  area 
to  that  of  the  cross-sections  is  given  by 

— V"  - 1.5  (1) 

2nr  / 3 

Since  the  above  result  is  Independent  of  particle  radius, 
the  exposed  surface  is  about  1.5  times  greater  than  that  of  a 
planar  surface.  Provisions  have  been  made  in  the  coda  for  ac- 
commodating such  a correction. 


TYPICAL  CRACK  THICKNESS  HISTORIES 


Limitations  of  One -Dimensional  Heat  Transfer  Analysis 

The  presence  of  surface  irregularities  raises  a question  as 
to  the  accuracy  of  treating  heat-transfer  into  the  propellant 
as  one-dimensional.  Of  particular  concern  is  whether  or  not  the 
irregularities  are  commensurate  with  the  depths  to  which  appre- 
ciable heat  is  transferred  into  the  propellant  and  the  depths 
of  material  evolved  during  critical  periods  of  burning. 

Depths  of  appreciable  heat  penetration  will  be  greatest  Just 
prior  to  ignition  and  can  be  crudely  estimated  using  the  expres- 
sion /at,  where  a represents  the  thermal  diffusivity  of  the  pro- 
pellant and  t represents  the  period  over  which  the  heat  is  being 
transferred.  Since  a is  in  the  vicinity  of  0.002  cm^/sec  for 
non-metalllc  ingredients,  and  the  heating  periods  can  be  as  small 
as  0.25msec  (see  Section  3),  it  is  concluded  that  the  depths  of 
appreciable  heat  penetration  are  of  the  order  of  7y  or  more  for 
the  clasB  of  problems  treated  in  this  study.  Depths  of  propellant 
evolved  during  burning  are  of  similar  magnitude. 

Unfortunately  the  depths  of  appreciable  heat  penetration 
are  not  always  small  compared  to  particle  dimensions,  Only  with 
large  particles  of  the  order  of  100m  or  greater  will  one -dimensional 
analysis  represent  an  adequate  approximation.  With  particles  of 
the  order  of  lOp  or  less,  one-dimensional  analysis  will  appreci- 
ably underpredict  the  velocities.  However,  until  the  significance 
of  such  errors  on  detonation  is  established,  it  is  premature  to 
embark  upon  a two-dimensional  analysis. 

2.3.2  Rates  of  Heat  Transfer  to  Propellant 

In  this  study  we  have  considered  the  convective  heating  of 
crack  surfaces  caused  by  the  flow  of  hot  high-velocity  high- 
pressure  gas.  Also  allowance  was  made  for  internal  heating  as- 
suming it  remains  constant  per  unit  mass  of  propellant  evolved. 


Because  of  uncertainties  in  the  thermal  gradients  within  the  gas, 
heat  feedback  from  burning  propellant  was  not  considered  other 
than  through  the  effect  of  burning  upon  the  bulk  temperatures, 
pressures  and  velocities  of  the  gases. 

Heat  fluxes  q are  evaluated  as  follows  s 

q - h(Tg(t)-T(t))  + Qp  p V (2) 

where  h - convective  heat-transfer  coefficient 

(see  Eq.  3) 

Tg(t) ,T(t)  - gas  and  propellant  temperatures  at  time  t 

Q - net  internal  heating  (reactions  plus  phase 
p changes)  per  unit  mass  of  evolved  propellant 

V,p  » propellant  velocity  and  density 

In  that  the  thickness  of  the  layer  in  which  internal  heating 
occurs  Is  very  small  (1)  compared  to  the  depth  of  heat  penetra- 
tion, the  net  internal  heating  is  assumed  to  occur  at  the  pro- 
pellant's surface. 

Net  Internal  Heating 

Indirect  determinations  (1)  of  the  magnitude  of  the  net  in- 
ternal heating  suggest  it  can  be  appreciable  for  some  secondary 
high  explosive  materials  during  steady  burning.  The  above  study 
suggests  values  as  high  as  130  cal/gm  for  HMX  during  low-pressure 
burning.  Values  estimated  by  other  investigators  for  other  pro- 
pellants are  of  similar  magnitude  (11,28,31).  Rates  of  pressure 
rise  have  important  implications  in  this  regard  through  their 
effect  upon  the  residence  times  of  reacting  gases. 

Convective  Heat-Transfer 

Convective  heat-transfer  coefficients  are  known  for  estab- 
lished flow  conditions  Over  given  surfaces.  Cracks,  of  course, 
present  highly  irregular  surfaces  to  the  gas  flow.  The  result 
is  large  uncertainties  in  the  flow  conditions  which  in  turn  af- 
fect the  heat- transfer  coefficient  h.  These  uncertainties  would 
be  greatest  during  the  final  stages  of  crack  collapse  wherein 


the  problem  becomes  most  complex  and  is  best  dealt  with  by 
experimental  means. 

As  a starting  point,  we  shall  neglect  local  variations  of 
the  heat-transfer  coefficient,  and  use  the  conventional  formu- 
lation for  the  heat  transfer  coefficient  h for  turbulent  flow 
given  by  McAdams  (34)  and  presented  below. 


K C Up  C 

h - ff  £ <-J4p£> 


0.8 


w 


(3) 


where  Kg  * thermal  conductivity  of  gas 

C„  - crack  thickness 
w 

U ■ gas  velocity 
Pg  - gas  density 

Cg^  - specific  heat  of  gas  at  constant  pressure 

C ■ dimensionless  constant  generally  taken  as  0.025 

The  constant  C is  inputted  into  the  code  to  allow  for  future  sen- 
sitivity studies  to  establish  the  importance  of  flux  uncertainties. 


2.3.3  Burn  Model 

An  important  consideration  in  modeling  burning  is  that  one 
must  account  for  highly  transient  temperatures  within  the  propel- 
lant involving  orders  of  magnitude  differences  of  time  for  igni- 
tion and  accelerated  burning.  To  expedite  the  calculations,  it 
was  necessary  to  develop  a novel  technique  for  predicting  the 
transient  temperatures  and  burning  rates.  This  involved  modify- 
ing an  existing  solution  of  the  temperatures  produced  in  a semi- 
infinite  body  by  a constant  flux  to  account  for  a moving  boundary 
subject  to  transient  fluxes.  Rates  at  which  gases  evolve  from 
the  propellant  were  made  by  applying  a first-order  Arrenhius  re- 
lationship to  the  propellant  temperatures  as  a function  of  depth. 
The  result  is  an  efficient  means  for  approximating  the  transient 
burning  rates. 


Two  versions  of  this  routine  were  compiled  --  one  for  exam- 
ining the  consequences  of  arbitrarily  chosen  time-dependent  fluxes 
upon  the  burning  rates,  and  one  for  determining  regression  veloc- 
ities at  various  locations  within  a crack. 

Means  for  Predicting  Transient  Propellant  Temperatures 

Predictions  of  transient  regression  rates  require  a knowledge 
of  the  transient  temperatures  produced  within  the  propellant.  Be- 
cause of  rapid  changes  in  the  heating,  quasi-steady  state  analysis 
will  not  suffice.  For  this  purpose  we  have  developed  an  approxi- 
mate solution  of  the  problem  involving  one-dimensional  heat  trans- 
fer into  the  propellant  with  a regressing  surface.  Means  for 
predicting  velocities  of  the  regressing  surface  are  discussed  in 
Appendix  D.  In  this  section  we  shall  brief ly  discuss  means  for 
predicting  the  transient  temperature  rises  and  refer  the  reader 
to  Appendix  D for  further  details, 

The  most  common  method  used  to  treat  thermal  problems  involv- 
ing regressing  surfaces  employs  finite-difference  techniques . 
Unfortunately  such  techniques  require  an  unacceptably  large  num- 
ber of  special  and  temporal  elements  to  follow  highly  dynamic 
burning  rates.  Moreover  no  closed-form  solution  of  the  problem 
exists . 

Hence  it  was  decided  to  develop  an  alternative  method  for 
predicting  the  propellant  temperatures.  This  method  is  predicated 
upon  constant  thermal  properties,  and  starts  with  the  solution  of 
the  problem  of  a semi-infinite  slab  exposed  to  a constant  heat 
flux.  Modifications  are  then  made  to  treat  the  problem  of  a 
regressing  surface  exposed  to  time -dependent  fluxes.  In  doing 
so  the  coordinate  system  is  fixed, and  the  regressing  surface 
moves  with  respect  to  the  fixed  coordinates.  Incremental  time 
intervals  are  used  during  each  of  which  the  incident  fluxes  are 
considered  constant. 

Essential  to  thi*  method  is  the  use  of  Duhamel's  principle 
to  determine  the  conductive  fluxes  produced  by  each  Incident  flux 


at  particular  depths  and  times  within  the  semi- infinite  body. 
These  depths  represent  given  depths  of  the  boundary  during  each  - 
of  the  time  intervals.  Negative  fluxes  are  introduced  to  elimin- 
ate conductive  fluxes  negated  by  the  loss  of  material.  By  doing 
so,  one  arrives  at  a number  of  constant  fluxes  at  the  given 
depths  with  which  to  approximate  the  heating  of  the  underlying 
propellant  over  time.  These  constant  fluxes  are  maintained  only 
over  the  time  interval  associated  with  the  particular  mean  depth. 

The  result  is  a number  of  simple  independent  heat-conduction 
problems  for  each  of  the  constant  fluxes.  Solutions  of  these 
problems  are  acquired  using  Duhamel's  principle  to  account  for 
differences  in  the  depths  and  times  at  which  the  various  fluxes 
are  applied.  Summing  the  various  temperature  contributions  one 
arrives  at  an  approximation  of  the  temperature  rise  of  the  bound- 
ary at  appropriate  times  caused  by  heat  conduction.  Here  it  is 
only  necessary  to  compute  surface  temperatures  to  predict  regres- 
sion rates,  even  though  the  method  could  have  just  as  easily  pre- 
dicted indepth  temperatures. 

Validation  of  Means  for  Predicting  Propellant  Temperatures 

To  check  the  above  procedure  for  predicting  transient  tem- 
peratures, we  have  applied  it  to  a problem  involving  a moving 
boundary  having  an  analytical  solution.  This  problem  is  presented 
below: 


£t.V8T.13T.0  x > 0 t >0 

^7  a a st  * * v c > u 

(A) 

T ■ TQ(  x>0,  t - 0 

(5) 

|~  - Ht,  x - 0 

(6) 

V - VoA  t > 0 

(7) 

where  T - temperature  at  time  t at  depth  x 
beneath  moving  boundary  at  x « 0 

a • thermal  diff us ivity 


u * heat  transfer  coefficient  divided 
by  thermal  conductivity 

V0  ■ velocity  of  moving  boundary  at  x ■ 0 


The  above  equations  describe  the  consequence  of  conduction  within 
a semi-infinite  body  at  a uniform  initial  temperature  with  a moving 
boundary  that  loses  heat  convectlvely  to  a zero  temperature  environ- 
ment. The  body  is  considered  to  move  at  a velocity  VQ  with  respect 
to  a fixed  coordinate  system  in  such  a way  that  the  moving  boundary 
is  always  located  at  x ■ 0.  Solution  of  the  above  problem  is  pre- 
sented on  page  389  of  Ref.  35  as 


T - 
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T -0.5  T (erfc  — 

° 0 2 /at 
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„c  x+V  t 

■ - exp  V x/a  erfc  — 
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x+(2aFi-V0)t 


> 


exp  (Fx-Kt  (V  -aH) ) erfc  

2<afi-Vrt>  ° 2 /at 


(8) 


Constants  used  were: 

2 

a - .0018  cm  /sec 
FT  - 277/cm 

TQ  - 100°C 

V “ 1,  3 and  5 cm/sec 

A comparison  of  the  temperatures  predicted  by  the  two  methods  is 
presented  in  Figs.  11,  12  and  13  for  the  velocities  cited  above. 

It  may  be  observed  that  the  temperatures  predicted  by  the  approxi- 
mate solution  of  Appendix  D are  in  good  agreement  with  those  pre- 
dicted by  the  exact  solution  presented  by  Eq.  8 for  all  three  of 
the  velocities  considered.  Thirty  time  Intervals  were  used  for 
each  of  the  approximate  solutions  — each  of  which  was  intention- 
ally kept  large  to  severely  test  the  approximate  method. 

The  reason  for  the  discrepancies  being  larger  at  the  moving 
boundary  is  due  to  the  fact  that  errors  in  the  stepwise  approximations 
of  the  incident  fluxes  having  their  greatest  impact  upon  the  sur- 
face temperatures.  Errors  in  the  changes  of  surface  temperature 
range  from  1.4  to  4.0  percent  with  the  highest  values  occurring 
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(10  time  intervals) 
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(20  time  intervals) 
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(30  time  intervals) 
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Figure  U TEMPERATURES  COMPUTED  BY  APPROXIMATE  AND  EXACT 
SOLUTIONS  FOR  A VELOCITY  OF  1 CM/SEC 


Temperature 


Elapsed  Time  N 0.10  msec 
(10  time  intervals) 

Elapsed  Time  ■ 0,30  msec 
(20  time  intervals) 

Elapsed  Time  ■ 0,54  msec 
(30  time  intervals) 


Note : Exact  solutions  given 
by  curves , approximat 
solution  by  dots 
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Figure  12  TEMPERATURES  COMPUTED  BY  APPROXIMATE  AND  EXACT 
SOLUTIONS  FOR  A VELOCITY  OF  3 CM/SEC 
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Figurt  13  TEMPERATURES  COftPUTED  BY  APPROXIMATE  AND  EXACT 
SOLUTIONS  FOR  A VELOCITY  OF  5 CM/ SEC 
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at  the  early  times  wherein  only  five  time  intervals  are  used. 
Thereafter  the  errors  decrease  progressively  with  time  to  values 
ranging  from  1 to  2 percent. 

Beneath  the  boundary  the  temperatures  are  in  good  agreement. 

In  this  regard  errors  associated  with  the  losa  of  heat  content 
are  less  than  one  percent  for  each  of  the  three  cases  presented 
in  Figs.  11,  12  and  13. 

Tabular  results  associated  with  the  three  figures  are  pre- 
sented in  Appendix  D.  Included  in  this  appendix  is  the  consequence 
of  varying  the  time  intervals  and  depths  at  which  the  fluxes  are 
applied  within  the  boundary  displacements. 

Several  observations  may  be  gained  by  examining  the  tabular 
results.  The  first  is  that  the  errors  decrease  in  a less  than 
linear  fashion  with  the  time  Intervals.  The  second  is  that  the 
errors  are  minimal  when  the  fluxes  are  applied  at  depths  corres- 
ponding to  approximately  three  tenths  of  the  displacements  of 
the  boundary.  Finally,  the  errors  are  relatively  insensitive  to 
the  velocity,  and  as  with  any  numerical  method  can  be  held  to  any 
desired  value  by  appropriate  selection  of  the  time  intervals. 

Means  for  determining  the  time  intervals  are  presented  in  Appen- 
dix D in  terms  of  the  velocity  at  which  the  propellants'  surfaces 
regresses;  accuracy  is  maintained  by  decreasing  the  time  Intervals 
as  the  velocity  increased. 

Means  for  Predicting  Regression  Ratea 

There  are  two  ways  that  propellants  can  regress  into  the  gas 
stream,  namely  by  the  discharge  of  molten  material  and/or  by  the 
evolution  of  gas.  To  appreciate  the  significance  of  each  phenom- 
enon, provisions  have  been  made  in  the  codes  to  treat  each  of 
these  ideal  situations  separately,  and  will  be  discussed  in  this 
section, 

Rates  at  which  gases  evolve  from  propellants  depend  primarily 
upon  the  temperature  distribution  within  the  propellant,  and  to  a 


lesser  extent  upon  the  pressure  to  which  the  propellant  is  sub- 
jected. Here  we  shall  neglect  the  effects  of  pressure  apart  from 
their  effect  in  increasing  the  flux.  Selection  of  the  most  appro- 
priate means  for  predicting  regression  rates  depends  upon  whether 
the  mass  of  critically  heated  propellant  remains  constant  or 
changes  in  an  inverse  fashion  with  the  flux.  For  situations  in 
which  the  mass  of  critically  heated  propellant  remains  essentially 
constant  the  following  expression  provides  a valid  approximation 
and  is  commonly  used: 

V - a exp-(E/(T(0)+To))  (9) 

where  a - constant 

E - activation  energy  divided  by  gas  constant 
T(0)  ■ temperature  rise  of  surface  of  propellant 

Tq  - initial  absolute  temperature  of  propellant 

For  situations  in  which  the  masa  of  critically  heated  pro- 
pellant decreases  in  an  inverse  fashion  with  flux,  it  is  more 
appropriate  to  integrate  the  Arrenhius  relationship  over  depth  x 
to  arrive  at  the  following  expression. 

V - Z J exp-(E/(T(x)+T0))dx  (10) 

*o 

where  Z - frequency  constant 

Means  for  approximating  the  above  integral  are  described  in  Ap- 
pendix A.  In  Appendix  A we  have  taken  advantage  of  the  fact  that 
most  of  the  gas  evolves  from  shallow  depths  compared  to  the  depths 
of  heated. propellant.  Temperature  rises  near  the  surface  are  as- 
sumed to  decrease  exponentially  with  depth  as  follows 

T(x)  - T (0)  exp-(Cx)  (11) 

Using  the  above  temperature  profile  in  conjunction  with  the 
thermal  gradient  at  the  surface,  the  integral  of  Eq.  10  is  approx- 
imated as  follows; 

V *KZ  exp-(E/(T(0)+To))(T(0)+To)2/(qff)  (12) 

where  K - thermal  conductivity  of  propellant 

q - effective  thermal  flux  (incident  plus  internal) 
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The  advantage  of  Eq.  12  oven  Eq.  10  is  in  eliminating  the  need 
to  calculate  temperatures  below  the  surface. 

For  situations  in  which  the  surface  regresses  due  to  the 
discharge  of  molten  propellant,  the  surface  temperature  is  main- 
tained at  the  melt  temperature  as  long  as  melting  occurB.  Pro- 
visions for  achieving  this  end  are  discussed  in  Appendix  A. 

In  both  of  these  highly  idealized  situations,  displacements 
of  the  moving  boundary  are  found  by  trial  and  error  methods . 

These  methods  are  discussed  in  Section  3 of  Appendix  A. 

2.3.4  Model  Results  for  Arbitrary  Thermal  Fluxes 

Here  we  shall  consider  the  effects  of  sudden  flux  increases 
upon  the  rate  of  surface  regression.  Figure  14  presents  the 
calculated  results  for  HMX  propellant,  while  Fig.  15  presents 
similar  results  for  a propellant  with  a lower  activation  energy  E. 

In  both  cases,  the  propellant  is  considered  to  be  initially 

o 

exposed  to  a constant  flux  of  500  cal/cm  -sec  which  is  maintained 
until  near  steady  velocities  are  achieved.  This  flux  is  similar 
in  value  to  convective  fixates  produced  by  hot  high-velocity  high- 
pressure  gases  in  cracks.  Then  the  flux  is  abruptly  increased 
by  an  arbitrary  factor  of  5 and  maintained  at  its  new  value  for 
all  future  times.  Values  of  the  constants  for  HMX  were  acquired 
from  Ref.  36  and  are  considered  as  independent  of  temperature. 

It  may  be  observed  that  the  consequence  of  the  sudden  flux 
increase  Is  a rapid  acceleration  of  the  velocity  to  high  values 
followed  by  a gradual  decrease  to  values  approaching  steady-state 
velocities.  The  overshoot  is  caused  by  the  sensible  heat  within 
the  propellant  being  greater  than  that  required  for  steady  burn- 
ing at  the  higher  flux.  Velocities  will  continue  to  exceed  their 
quasi-steady  state  counterparts  as  long  as  the  flux  continues  to 
rise.  If  the  flux  remains  constant,  the  excess  sensible  heat 
will  be  consumed  by  the  burning  and  the  velocity  will  decrease 
asymptotically  to  steady-state  values  as  shown  in  Figs.  14  and  15. 
If  the  flux  decrease's  appreciably,  the  velocities  will  drop  below 
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Figure  14  CALCULATED  EFFECT  OF  ABRUPT  RISE  OF  HEAT  FLUX 
UPON  REGRESSION  VELOCITY  OF  HMX 


Figure  15  CALCULATED  EFFECT  OF  ABRUPT  RISE  OF  HEAT 
UPON  REGRESSION  VELOCITY  OF  PROPELLANT  W 
E/R  * 15 ,Q00°K 


their  quasi-steady-state  values  and  could  extinguish  further 
burning.  This,  of  course,  is  the  reason  for  the  recovery  of  the 
explosive  fragments  shown  in  Fig.  1. 

It  may  be  observed  that  the  overshoots  are  most  pronounced 
with  the  higher  E/R  value.  For  example,  the  peak  velocity  of 
Fig.  14  is  about  3 times  greater  than  their  quasi-steady-state 
counterpart  while  Fig.  15  indicates  a factor  of  about  2. 

In  addition  to  dampening  any  velocity  increases,  lower  E/R 
values  also  delay  the  occurrence  of  the  peak  velocities.  This 
could  have  important  implications  in  the  problem  of  crack  burning 
since  the  delay  times  are  of  similar  magnitude  as  the  time  in- 
volved in  the  final  stages  of  crack  collapse  wherein  high  pressures 
are  expected. 

Each  of  the  above  effects  is  not  at  all  surprising  in  that 
propellants  with  higher  E/R  values  are  more  sensitive  to  tem- 
perature charges.  The  primary  factor  here  is  differences  in 
the  temperature  rises  necessary  to  achieve  specific  velocity 
changes.  These  temperature  rises  are  less  with  propellants  having 
higher  E/R  values. 

Another  possibly  important  factor  in  accentuating  the  veloc- 
ities is  internal  heating.  Exothermic  heating  over  and  above  that 
needed  for  sublimation  will  accentuate  each  of  the  above  effects  -- 
namely  more  rapid  and  pronounced  velocity  increases.  However, 
until  more  is  known  regarding  the  Internal  heating  produced  during 
extremely  transient  high-pressure  burning,  one  cannot  quantify 
this  effect. 


3.  DYNAMIC  MODEL  OF  CRACK  BURNING 


This  model  assesses  the  transient  burning  rates  and  pressures 
caused  by  the  sudden  exposure  of  a crack  (within  a monopropellant) 
to  a cavity  of  hot  gases  held  at  a constant  pressure  and  temperature. 
The  initial  pressure  and  temperature  of  the  gas  within  this  crack 
is  considered  much  lower  than  that  of  the  cavity  in  order  to  accen- 
tuate the  burning  and  crack  deformations,  and  in  turn  the  possi- 
bilities for  pronounced  pressure  buildup  within  the  crack. 

Major  assumptions  used  to  expedite  development  of  the  model 

were 

e One -dimensional  gas  flow  with  friction 

e Perfect  gas 

e Instantaneous  reaction  of  evolved  propellant 

e Propellant  behaves  as  a simple  isotropic  visco- 
elastic material 

e Fully  established  turbulent  gas  flow 

e One-dimensional  heating  of  propellant  surfaces 

e Constant  gas  and  propellant  properties 
3 . 1 Computer  Code 

This  code  has  been  designed  primarily  as  a research  tool  to 
assist  in  Identifying  the  essential  features  of  DDT.  As  a con- 
sequence it  is  subject  to  periodic  revisions  as  key  aspects  of 
the  problem  are  identified.  For  this  reason  we  shall  limit  this 
discussion  to  the  principal  features  of  the  code.  Specific  tech- 
niques used  to  treat  various  aspects  of  the  problem  are  discussed 
elsewhere  in  this  report,  Functions  performed  by  various  portions 
of  the  code  are  described  by  Fig,.  16. 

The  main  program  and  each  of  the  subroutines  are  used  once 
during  each  time  Interval.  These  Intervals  are  computed  in  the 
main  program  using  criteria  involving  the 


MAIN  PROGRAM 

Computes : 

• Time  intervals 

• Propellant  temperattures/velocities 
associated  with  each  crack  element 

• Rates  of  energy/mass  flows  to  gas 
stream 

• Changes  in  crack  thicknesses  due 
to  burning 

• Arrays  needed  to  consider  crack 
propagation 


Figure  16  PRINCIPAL  FEATURES  OF  DYNAMIC  MODEL 


• Propellant  regression  rates  along  with  thermal 
fluxes  (see  Appendix  A) 

• Gas  and  sound  velocities  (see  Appendix  B) 

e Crack  thicknesses  and  deformation  velocities 
(see  Appendix  C) 

These  criteria  are  applied  to  each  of  the  crack  elements  to  find 
the  smallest  time  interval  required.  In  this  regard,  the  UBer 
has  the  option  of  using  somewhat  larger  time  intervals  than  needed 
by  the  gas  dynamic  calculations  whenever  the  consequence  upon  the 
propellant  is  small.  Time  intervals  decrease  as  the  velocities 
associated  with  the  gas,  propellant,  and  crack  deformations  increase. 

Following  the  selection  of  the  time  interval,  subroutine 
Gasdyn  is  called  to  compute  the  gas  pressures,  velocities  and 
temperatures  associated  with  each  crack  element.  This  subroutine 
also  calculates  the  mechanical  deformations  of  the  crack. 

Subroutine  Gheat  then  uses  the  predicted  gas  pressures,  tem- 
peratures  and  velocities  to  compute  the  convective  fluxes  entering 
the  propellant  at  each  crack  element.  These  fluxes  are  then  added 
to  the  net  internal  heating  to  arrive  at  the  total  flux.  Finally 
the  main  program  is  reentered  to  compute  the  propellant  tempera- 
tures/velocities . Before  commencing  the  next  time  interval,  the 
rates  of  energy  and  mass  discharge  from  the  propellant  are  cal- 
culated along  with  changes  in  crack  thickness  due  to  propellant 
regression  for  subsequent  use  by  subroutine  Gasdyn. 

As  mentioned  previously,  provisions  have  been  made  in  the 
code  for  investigating  the  consequence  of  various  aspects  of  the 
problem.  These  include 

e Variable  crack  shapes,  lengths,  thicknesses  and 
end  conditions 

e Crack  propagation  (at  present  induced  arbitrarily) 

e Double  or  single  propellant  surfaces  to  distinguish 
between  propellant  crack  from  motor  case  debonds 


• Choice  of  one  of  two  idealized  situations  in  which 
molten  propellant 

- remains  in  place 

- is  immediately  swept  into  gas  stream 
upon  formation 

3.2  Results  from  Dynamic  Model 

In  this  section  we  shall  present  some  of  the  more  pertinent 
findings  obtained  from  dynamic  model.  In  view  of  various  simpll- 
fications  and  uncertainties  in  the  data,  the  results  are  of  a 
preliminary  nature  and  are  primarily  intended  to  identify  essen- 
tial features  of  the  problem  requiring  more  detailed  considera- 
tion and  study.  In  each  of  the  cases  treated,  the  propellant 
surfaces  are  assumed  to  regress  by  the  evolution  of  gas . 

Tables  1 and  2 indicate  the  effects  of  cavity  pressure,  and 
crack  thicknesses  and  end  condition  upon  the  spread  of  burning 
and  upon  the  mechanical  deformation  of  cracks.  In  each  of  these 
tables  the  downstream  end  of  open  cracks  is  considered  connected 
to  a constant-pressure  reservoir  held  at  ambient  pressure.  Burn- 
ing is  arbitrarily  defined  as  occurring  when  the  regression  veloc- 
ities exceed  0.1  cm/sec. 

Tables  1 and  2 Indicate  that  the  effect  of  cavity  pressure 
upon  the  spread  of  burning  is  much  more  than  linear.  Also  it 
may  be  observed  that  high  upstream  pressures  tend  to  open  the  up- 
stream end  of  the  crack  and  close  the  downstream  end  of  the  crack. 
The  degree  to  which  various  portions  of  a crack  open  or  close  de- 
pends upon  the  gas  pressures  within  the  crack  and  the  intensity 
of  the  reflected  stress  wave.  In  this  regard,  stiff  motor  cases 
in  close  proximity  to  the  crack  are  most  conducive  to  crack  closure. 

Over  the  range  of  parameters  considered,  Tables  1 and  2 in- 
dicate that  crack  thickness  affects  the  spread  of  burning  less 
than  the  cavity  pressure.  This  may  not  be  true  at  lower  gas 
pressures  and  crack  thicknesses.  Also  it  may  be  observed  that 
with  long  cracks  of  the  order  of  feet  in  length  that  the  down- 
stream end  condition  has  a negligible  effect  upon  the  initial 
spread  of  burning.  In  this  regard  experimental  studies  (14) 


Table  1 


SPREAD  OF  BURNING  FROM  CAVITY  INTO  CRACKS 
AS  FUNCTION  OF  CAVITY  PRESSURE* 


Cavity 

Pressure, 

Bars 

Depths  of  Burning  as  Function  of  Time, 
0.3  ms  0.4  ms  0.5  ms  0.6  ms 

cm 

0.7  ms 

50 

0 

(.33-. 30) 

0 

(.34-. 29) 

0.5 

(.34-. 29) 

1.0 

(.34-. 28) 

“To" 

(.34-. 27) 

100 

8.0 

(.36-. 30) 

16.0 

(.37-. 29) 

20.5 

(.38-. 28) 

22.0 

(.38-. 26) 

26.5 

(.39-. 22) 

200 

28.0 

(.42-. 20) 

34.0 

(.45-. 28) 

38.0 

(.47-. 26) 

40.0 

(.49-. 22) 

40.5 

(.50-. 18) 

* Cavity  gasat  2300°K;  uniform  cracks  0.30  cm  wide,  60  cm  long, 
open  end 

**  Numbers  in  parenthesis  refer  to  crack  thicknesses  in  cm  at 
front  and  rear  of  crack  as  influenced  by  gas  pressures  and 
reflected  stress  waves  (case  stiffness  and  distance  of  motor 
case  assumed  as  5,000  Bars /cm  and  5.  cm,  respectively) 


Table  2 

SPREAD  OF  BURNING  FROM  CAVITY  INTO  CRACKS 
AS  FUNCTION  OF  INITIAL  CRACK  THICKNESS* 


Initial 
Crack 
Width,  cm 

Depths  of  Burning  as  Function  of  Time 

cm 

0.3  ms 

0,4  ms 

0.5  ms 

0 .6  ms 

0.7  ms 

0.05 

7.0 

8.5 

9.5 

10.5 

11.5 

(open) 

(.11-. 05) 

(.12-. 04) 

cn 

o 

1 

<r 

rl 

(.16-. 01) 

(.17-0) 

0.12 

7,5 

12.0 

15.0 

16.0 

17.0 

(open) 

(.36-. 30) 

(.37-. 29) 

(.38-. 28) 

(.38-. 26) 

(.39-. 22) 

0.30 

8.0 

16.0 

20.5 

22.0 

26.5 

(open) 

o* 

1 

u 

o 

(.37-. 29) 

(.38-. 28) 

(.38-. 26) 

(.39-. 22) 

0.30 

8.0 

16.0 

20.5 

21.5 

25.5 

(closed) 

(.36-. 29) 

(.37-. 29) 

(.38-. 29) 

(.38-. 30) 

(.38-. 32) 

* Caylty  gas  at  200  Bars,  2300  K;  uniform  cracks,  60  cm  long 
**  Numbers  in  parenthesis  refer  to  crack  thicknesses  in  cm  ac 
front  and  rear  of  crack  a*  influenced  by  gas  pressures  and 
reflected  stress  waves  (case  stiffness  and  distance  of  motor 
case  assumed  as  5,000  Bars /cm  and  5.  cm,  respectively) 
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with  short  cracks  of  the  order  of  inches  in  length  show  that  the 
end  condition  has  a very  important  effect  upon  the  spread  of  burn 
ing  under  essentially  uniform  pressures. 

Finally  it  should  be  observed  that  the  downstream  end  of 
closed-ended  cracks  tend  to  contract  less  than  cracks  connected 
to  a low-pressure  reservoir.  In  none  of  the  above  cases  did  the 
gas  pressures  within  the  crack  rise  appreciably  above  the  cavity 
pressure. 

In  order  to  produce  pronounced  pressure  rises  within  the 
crack  we  chose  a long  crack  with  an  inverted  taper  to  yield  more 
uniform  crack  thicknesses  following  contraction.  An  inverted 
taper  was  chosen  in  that  cracks  tend  to  open  at  the  upstream  end 
and  close  at  the  downstream  end.  The  upstream  end  of  the  crack 
was  considered  connected  to  a cavity  held  at  a pressure  of  200 
Bars  while  the  downstream  end  was  considered  closed.  To  accen- 
tuate the  reflected  stress  wave  we  chose  a very  stiff  motor  case 
of  10,000  Bars/cm  located  at  a distance  of  5 cm  from  the  crack. 

Crack  deformations,  and  gas  pressure  and  velocities  are 
shown  in  Fig.  17.  The  crack  deformations  are  scaled  independ- 
ently of  the  crack  length.  Extent  of  burning  is  indicated  by 
the  dashed  lines  using  an  arbitrary  burn  criterion  of  0.1  cm/sec. 
Gas  velocities  and  pressures  within  the  crack  ara  presented  along 
the  length  of  the  crack.  Pressure  and  velocity  pre'dlctlono  are 
not  presented  at  0,73  ms  due  to  unaccountable  energy  losses  from 
the  gas.  It  is  believed  that  tha  computational  error  is  due  to 
order  of  magnitude  variations  in  tht  crack  thloknasses  during 
the  final  stages  of  crack  collapse.  Throughout  the  run,  mass  was 
conserved.  In  spite  of  large  energy  losses  very  rapid  and  pro- 
nounced  pressure  rises  were  predicted  within  the  crack  at  0.73  ms. 
However,  the  magnitude  of  them  pressures  cannot  be  defined  until 
the  ooda  Is  revised  to  treat  pronounced  variations  in  tha  crack 
thicknesses. 

From  Fig.  17  it  may  bs  seen  that  the  gas  prer suras  are 
generally  low  up  to  0,50  tbs  except  at  the  closed  end  of  tha  crack. 


Cavity;  20Q  Bara,  23Q0°K 

Crack  45  cm  long,  0.20-0.29  cm  thick 


20Q  Bara 


200  Bars 


0 cm/sec 

“ 

1 Bar 

70,700  cm/sec \ 

\ 12Q 

0 

1 

74, 100  5 

\ 100 

0 

1 

81,300 

77 

0 

1 

95,700 

59 

o 

1 

112,000 

48 

0 

1 

126,000 

41 

0 

1 

137,000 

37 

0 

1 

144,000 

• 

34 

0 

1 

144 , 000 

32 

0 

1 

119,000  i 

39 

Time-0  Time-,25  ms 


200  Bars 


200  Bara 


Time-. 7 3 ms 


Figure  17  SPREAD  OF  BURNING t CRACK  DEFORMATIONS { 
AND  OAS  PRESSURES  AND  VELOCITIES 
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The  elevated  pressure  is  caused  by  the  reflected  shock.  Up  to 
this  time  high  velocity  gas  flow  is  the  prime  cause  of  pressure 
changes  rather  than  burning.  The  reverse  would  be  expected  if 
the  crack  contracts  sufficiently.  Important  considerations  in 
this  regard  are  the  extent  to  which  a crack  is  able  to  contract 
in  view  of  its  highly  Irregular  surfaces , and  the  consequences 
of  possible  impact.  In  the  example  cited  by  Fig,  17,  the  minimum 
crack  thickness  was  achieved  at  0.73  ms  and  subsequently  increased. 
However,  in  view  of  the  computational  error  cited  earlier,  it  was 
not  possible  to  quantify  how  close  the  propellant  surfaces  came 
to  impacting. 

Depths  of  propellant  evolved  were  quite  small  --  ranging 
from  a few  microns  to  a few  tens  of  microns.  The  greatest  loss 
occurred  at  the  upstream  end  of  the  crack  due  to  higher  heat- 
transfer  rates  and  longer  burning  times.  At  the  region  of  mini- 
mum crack  thickness  at  0.73  ms,  only  2 microns  of  propellant 
were  evolved.  This  suggests  that  crack  surfaces  may  not  be  ap- 
preciably altered  by  the  burning  prior  to  any  impact.  Leaving 
much  of  the  more  sensitive  propellant  components  at  the  surface 
could  have  important  effects  upon  the  consequences  of  crack 
closure/impact.  This  subject  will  be  returned  to  in  the  next 
section. 
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4.  SUMMARY,  CONCLUSIONS  AND  FUTURE  NEEDS 

In  this  section  we  shall  summarize  conclusions  derived  from 
this  study.  Following  this  we  shall  present  our  views  regarding 
future  studies  to  resolve  this  most  complex  problem. 

4.1  Summary  and  Conclusions 

Several  phenomena  have  been  identified  as  playing  an  impor- 
tant role  in  causing  rapid  pressure  rises  within  cracks.  The 
first  is  pressure-induced  deformations  of  cracks. 

During  periods  of  gas  filling,  cracks  will  expand  in  a non- 
uniform  fashion  as  high-pressure  gas  flows  into  the  crack.  Crack 
surfaces  can  displace  outward  at  rates  as  high  as  100  to  200  cm/sec. 
Since  the  period  of  gas  filling  is  of  the  order  of  milliseconds 
(depending  upon  crack  dimensions) , cracks  can  expand  by  as  much 
as  tenths  of  a centimeter. 


Cracks  will  continue  to  expand  at  least  until  arrival  of 
stress  waves.  These  stress  waves  could  arise  from  other  cracks 
and/or  by  reflection  from  tha  motor  case.  In  this  study  con- 
sideration was  given  only  to  reflected  stress  waves.  With  mul- 
tiple cracks,  the  problem  becomes  highly  statistical.  Quanti- 
fication of  the  conditions  (cracks,  motor  case,  and  propellant 
composition)  necessary  to  produce  rapid  closure  of  cracks 
has  yet  to  be  accomplished. 

Another  possibly  important  factor  in  causing  rapid  pressure 
rises  is  that  burning  velocities  can  achieve  values  at  least  2 
to  3 times  larger  than  their  quasl-steady-state  counterparts  as 
a consequence  of  rapid  increases  of  heating.  In  this  regard, 
increased  gas  pressures,  temperatures  and  velocities  produced 
during  the  final  stages  of  crack  collapse  could  have  important 
implications  in  accelerating  gas  burning. 

At  least  two  aspects  of  propellants  can  play  an  important 
role  in  accentuating  the  propellant  burning  rates.  The  first 
is  their  activation  energy.  The  Higher  the  activation  energy, 


46 


the  less  the  temperature  rise  needs  to  be  to  increase  the  velocity 
by  a given  amount.  Thus  propellants  with  higher  activation  ener- 
gies will  exhibit  more  rapid  and  pronounced  velocity  responses 
to  changing  thermal  fluxes.  In  this  regard  HMX  is  particularly 
hazardous,  in  that  its  activation  energy  is  almost  twice  that 
of  conventional  propellant  ingredients. 

The  second  factor  is  internal  heating.  During  the  course 
of  this  study  we  have  found  that  positive  net  internal  heating 
will  cause  greater  overshoots  of  the  velocities  than  those  in- 
dicated previously  for  zero  net  internal  heating  (see  Figs.  14 
and  15) . The  problem  here  is  the  lack  of  information  describing 
internal  heating  during  highly  dynamic  burning.  Using  steady- 
state  values,  while  questionable,  is  the  best  alternative  at 
the  present  time. 

During  the  final  stages  of  crack  collapse  the  problem  be- 
comes more  and  more  conjectural  involving  large  uncertainties 
in  the 

• Consequence  of  possible  impact 

• Gas  flow 

• Amounts  of  heat  feedback  from  burning  propellant 

e Internal  heating 

• Reaction  kinetics 

At  some  point  during  crack  collapse,  analysis  will  not  suf- 
fice and  one  must  resort  to  experiments.  Model  predictions  will 
still  be  needed,  however,  to  establish  the  experimental  conditions. 
In  the  next  section  we  shall  briefly  discuss  the  type  of  experi- 
ments needed  to  assess  the  consequence  of  crack  collapse. 

4.2  Future  Needs 

Several  important  provisions  need  to  be  incorporated  into 
the  code.  These  are 

e Equation  of  state 

e Variable- length  crack  elements  to  expedite  calculations 


• Means  for  approximating  heat  feedback  from  burning 
propellant 

• Means  for  accomodating  highly  variable  crack 
thicknesses  in  the  gas  dynamic  calculations 

• Crack  propagation  criterion 

• Provisions  for  branching  cracks 

• Time-dependent  conditions  at  open  ends  of  cracks 

• Detonation  criterion 

Of  paramount  importance  is  knowledge  regarding  the  consequence  of 
impacting  burning  surface  at  elevated  pressures.  As  noted  earlier, 
such  knowledge  is  best  acquired  experimentally.  These  experiments 
should  involve  impacting  burning  propellant  surfaces  ar  various 
ambient  pressures,  In  such  experiments  care  needs  to  be  exercised 
to  achieve  realistic  closure  speeds,  as  well  as  crack  surfaces 
and  burning  conditions. 

The  ultimate  objective  of  such  experiments  should  be  the 
identification  of  the  threshold  conditions  (closure  speeds,  burn- 
ing conditions  and  pressures)  required  to  initiate  det. nation. 

Once  these  threshold  conditions  have  been  established,  the  code 
Bhall  be  in  a position  to  explore  various  means  by  which  the 
threshold  conditions  may  be  avoided  through  changes  in  high-energy 
propellants,  operational  pressures,  and  motor  design. 
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APPENDIX  A 


MEANS  FOR  PREDICTING  REGRESSION  VELOCITIES 
A. 1 INTRODUCTION 

During  burning,  propellants  can  evolve  material  into  the 
gas  stream  by  the  discharge  of  molten  material  and  gas.  These 
effects  arc  not  independent  of  each  other  in  that  molten  mate- 
rial from  one  of  the  propellant  constituents  can  flow  over  other 
constituents  and  affect  the  rates  at  which  they  evolve.  Here 
we  shall  discount  such  an  effect  and  concentrate  our  attention 
upon  each  of  the  two  ideal  situations  cited  above.  The  actual 
situations  undoubtedly  lie  somewhere  in  between. 

A.?.  CRITERIA  FOR  ASSESSING  VELOCITIES 

A. 2.1  Gas  Evolution 

A common  practice  is  to  predict  regression  velocities  V 
using  the  surface  temperature  T(0)  +TQ  as  follows 

V - a exp  - <Er/(T(0)+To)  <A-1) 

where  a - constant,  cm/sec 

E * activation  energy  divided  by  gas  constant 

Implicit  in  this  expression  is  that  there  is  a constant  mass  of 
material  at  temperature  T(0)  + TQ.  Here  we  refer  to  the  presence 
of  a thin  hot  melt/solid  layer  from  which  much  of  the  gas  origin- 
ates. Certainly  this  assumption  is  reasonable  for  steady  or 
near-steady  burning  over  a limited  range  of  velocities. 

In  order  to  use  Eq.  A-l  it  is  first  necessary  to  predict 
the  surface  temperature  T(0).  This  is  usually  accomplished  by 
assuming 

e all  heat  transfer  through  «ny  molten  layer  and 
underlying  propellant  is  conductive 

e all  internal  heating  and  latent-heat  absorption 
occurs  at  the  surface 


While  neither  of  these  assumptions  are  appealing  they  are  pres- 
ently necessary  in  lieu  of  adequate  information  and  understand- 
ing of  the  many  complex  processes  taking  place  near  the  surface. 

In  the  remainder  of  this  appendix  we  shall  present  an  altern- 
ative to  Eq.  A-l,  predicated  upon  changes  of  the  mass  from  which 
gas  evolves  during  highly  dynamic  burning.  Choice  of  most  appro- 
priate approach  depends  upon  how  the  mass  of  critically  heated 
propellant  changes  during  dynamic  heating  and  burning.  With  Com- 
position B,we  have  found  the  foam  mass  during  steady  burning  to 
increase  slightly  with  pressure  (1).  On  the  other  hand,  PBX  re- 
sulted in  decreases  in  mass  as  the  pressure  was  elevated.  Unfor- 
tunately such  steady-state  experiments  do  not  Indicate  how  the  mass 
will  change  under  dynamic  heating  conditions. 

For  this  case  we  shall  assume  that  the  melt  is  indispersed 
with  solid  particulates  in  such  a fashion  that  one  can  treat 
the  heat  transfer  as  conductive.  All  thermal  properties  are 
assumed  as  constant.  Integrating  the  Arrenhlus  relationship 
over  some  depth  6 containing  most  of  the  heat  yields 

V - Z f exp  (-E/(T(x)+T  ))dx  (A-2) 

where  V « rate  of  surface  regression 
T(x)  - temperature  rise  at  depth  x 

Z - frequency  factor 

£ - activation  energy  divided  by  universal  gas  constant 
Tq  ■ initial  temperature,  absolute 

Because  of  the  exponential  term  in  Eq.  A-2  the  Integration 
needs  to  be  performed  only  over  shallow  depths  of  propellant 
having  temperatures  in  the  vicinity  of  the  surface  temperature. 
Proceeding  upon  this  fact,  the  near-surface  temperatures  T(x) 
will  be  approximated  using  the  following  exponential  relationship. 

T(x)  - T(0)  exp-Cx  (A-3) 

where  C is  a constant  which  for  steady-state  burning  equals  the 
velocity  V divided  by  the  thermal  dlffusivity  a.  To  determine 

A3  , 


the  appropriate  expression  for  C for  unsteady  burning  we  shall 
multiply  Eq.  A-3  by  -K  and  differentiate  with  respect  to  x to 
arrive  at  an  effective  flux  q of 


q - -K  ^^|x.0  “ KT<°>C 


<A-4) 


where  K ■ thermal  conductivity.  Solving  Eq.  A-4  for  C yields 


c " rt for  <A_5) 

Equation  A-2  will  now  be  approximated  by  first  replacing 
the  variable  x by  the  variable  y given  below. 

y ■ T(0)  exp  <-Cx)  + TQ  (A-6) 


where  y2  - T<0)  + TQ 
yl  " T<6)  + TQ 

Expanding  the  integral  of  Eq.  A-7  by  parts  yields 


dy  - SS^SmpiL  |y2.2  fyZ^ln^X  dy  + 

J y-To  E(y-T0)  J E(y-T0) 

\ *1  (A-8) 

f1  »x£-(E/y)y2dy. 

JY1  *<y-Vz 

Here  the  two  integral*  on  the  right-hand  side  of  Eq.  A-8  are 
small  compared  to  the  integral  on  the  left-hand  side  and  tend 
to  cancel  each  other  out.  Neglecting  these  terms, and  the  negli- 
gibly small  contributicns  from  the  lower  limit  y^ , yields 

V.^  exp(-E/<T(0)+To).(T(0)+to)2/T(0)  (A-9) 

Using  the  appropriate  expression  for  C,  this  equation  can  be 
applied  to  either  steady  or  unsteady  burning.  For  unsteady  burn- 
ing, the  Velocity  from  Eqe.  A-5  and  A-9  Is 


(A- 10) 


V*  ~ exp(-E/(T(0)+T  ) ) • (T(0)+T  )2 
qE  oo 

This  equation  presumes  that  the  critical  mass  decreases  in- 
versely with  the  heat  flux.  Choice  of  Eq.  A-l  or  Eq.  A-l  depends 
upon  whether  or  not  this  assumption  is  more  valid  than  a constant 
mass.  Here  we  have  arbitrarily  chosen  to  use  Eq.  A-l  for  HMX. 
Future  analyses  needs  to  be  conducted  to  establish  the  consequence 
of  the  two  assumptions  upon  detonation.  If  found  important,  ex- 
perimental studies  should  be  conducted  to  establish  which  of  the 
two  approaches  is  more  valid  for  each  of  the  more  sensitive  pro- 
pellant ingredients. 

A. 2.2  Melt  Evolution 

To  treat  this  ideal  situation,  the  temperature  of  the  regres- 
sing surface  must  be  maintained  at  its  melt  temperature.  Incident 
fluxes  q 1 are  discounted  by  the  rate  of  heat  absorption  needed  to 
melt  the  propellant  at  the  computed  velocity.  The  result  is  repre 
sented  by  the  effective  flux  q.  Temperatures  are  calculated  as 
described  by  Appendix  D. 

In  the  next  section  we  shall  discuss  how  the  required  dis- 
placements are  calculated  for  each  time  Interval. 

A. 3 COMPUTATIONAL  SCHEMES  FOR  DETERMINING  DISPLACEMENTS  Axj 

Displacements  Ax  are  determined  for  each  time  interval  by 
a trial  and  error  method  using  the  criteria  discussed  in  Section 
2 of  this  appendix.  Here  the  displacement  or  velocity  ere  ap- 
proximated by  computing  limits  within  which  the  correct  value 
lies  — Ax^  for  the  case  of  gas  evolution,  of  for  the  case  of 
melt  evolution.  Theae  limits  are  progressively  narrowed  following 
each  trial.  Successive  estimates  are  set  close  to  the  geometric 
mean  for  the  case  of  gas  evolution,  or  equal  to  the  arithmetic 
mean  for  the  case  of  melt  loas. 

A. 3.1  Diaplacements  by  Gas  Evolution 

For  each  trial  the  temperature  at  x^  + Ax^  is  computed 
using  Eqs.  D-21  and  D-22  followed  by  a determination  of  the 
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velocity  using  Eq.  A-10.  The  estimated  Ax^  and  the  resultant 
value  for  Ax^  - (V^j+V^)  At^/2  are  then  used  to  narrow  the  ax 
limits.  If  Axi>Ax^,  the  lower  limit  is  set  equal  to  Ax^  while 
the  upper  limit  is  set  equal  to  0. SAXj+O . fiAx^ . 

If  Ix^  < Ax^  the  upper  lljnit  is  set  equal  to  Ax^  and  the 
lower  limit  is  set  equal  to  0.5Ax^+0.4Ax^. 

The  reason  for  setting  the  limits  as  described  Is  due  to 
the  fact  that  temperature  decreases  with  depth.  The  result  is 
non-linear  decreases  in  the  estimates  of  the  velocity  with  in- 
creases in  depth  such  that 

|Ax^-Ax|  > |Ax^-Ax|  (A-12) 

where  Ax^  and  Ax^  bracket  the  correct  value  Ax. 

If  the  limits  agree  to  within  5 percent  the  computations 
are  completed.  Otherwise,  a new  trial  Ax^  is  taken  cloBe  to 
the  geometric  mean  of  the  two  limits  and  the  process  repeated. 
Generally  2 to  7 trials  are  needed. 

A . 3 . 2 Displacements  by  Melt  Evolution 

This  situation  is  treated  in  a very  similar  fashion  as 
discussed  in  Section  3.1.  The  principal  differences  are  in 
using  velocity  limits,  and  in  using  arithmetic  means  to  upgrade 
the  velocity  estimates. 

For  each  trial  velocity,  the  trial  displacement  Is  substi- 
tuted into  Eqs . D-21  and  D-22  to  determine  the  surface  temperature 
Then  the  velocity  limits  are  narrowed  according  to  the  calculated 
temperature  T^  at  the  trial  depth  Ax^.  For  the  first  trial,  the 
upper  velocity  limit  is  set  equal  to  the  highest  value  possible 
namely 

Vplftl  <A'13) 

m 

where  q - effective  flux  allowing  for  melting  (See  Eq.  2) 
p - density 
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C - specific  heat 

Tm  * melt  temperature  above  initial  temperature 

Initially  the  lower  limit  is  set  equal  to  0.  A trial  velocity 
V is  then  taken  hAlfway  between  the  two  limits , and  used  to  deter- 
mine the  temperature  rise  corresponding  to  the  resultant  dis- 
placement V*  At.  Here  At  represents  the  time  Interval  previously 
selected. 


Limits  are  then  adjusted  according  to  whether  or  not  is 
above  or  below  T . Here  we  use  the  temperature  gradient  to  esti- 
mate the  displacement  and  velocity  changes  necessary  to  yield  Tm> 
To  prevent  overcorrections,  only  half  of  the  velocity  changes  are 
used.  The  result  is  the  new  limiting  velocity  7 presented  below 


7 


V + 


2qAt 


(A-14) 


If  T^  is  less  than  Tm,  the  upper  limit  is  set  equal  to  the 
minimum  of  its  previous  value  and  7.  If  T^  exceeds  Tm,  is 
set  equal  to  the  maximum  of  its  previous  value  and  7. 


A new  trial  vulue  la  then  computed  for  V based  upon  the 
arithmetic  mean  of  and  V^.  Once  these  limits  agree  to  within 
5 percent  the  computations  are  completed.  Generally  2 to  7 trials 
are  necessary  to  approximate  the  velocity  and  associated  displace- 
ment Ax. 


APPENDIX  B 
GAS  DYNAMIC  MODEL 


B.l  INTRODUCTION 

One  essential  aspect  of  the  initiation  and  burning  of  propel- 
lants in  cracks  is  the  transfer  of  the  gaseous  reaction  products 
(and  the  energy  and  momentum  associated  with  this  gas)  along  the 
crack.  Thus  the  establishment  of  a gas  dynamic  model  which  in- 
cludes a variety  of  flow  features  which  may  be  significant  is  a 
basic  requirement  of  this  study.  This  gas  dynamic  model,  when 
properly  coupled  with  the  heat  transfer,  ignition,  and  burn  models 
will  allow  for  the  appropriate  interaction  of  these  physical  events 
to  occur  and  permit  an  investigator  to  examine  the  severity  and 
sequence  of  the  events  which  occur  in  this  complex  problem  for  a 
variety  of  parameter  values  or  model  variations  or  assumptions. 

In  this  manner  an  insight  into  the  occurrence  of  the  accidental 
detonation  of  high  energy  propellants  can  be  developed  as  well 
as  providing  a rational  basis  for  the  design  of  experiments  with 
which  to  verify  or  evaluate  the  importance  of  some  parameters 
or  effects.  This  appendix  presents  the  development  of  the  gas 
dynamic  model  and  its  associated  numerical  representation  for 
computer  application. 

The  gas  dynamic  model  which  already  includes  a variety  of 
features  and  effects  is  still  considered  to  be  a preliminary 
version  subject  to  modifications.  These  modifications  will  be 
indicated  by  initial  results  and/or  the  limited  treatment  of 
some  aspects  of  the  Initial  model.  The  gas  dynamic  model  is 
based  upon  an  Eulerlan  representation  of  the  flow  region  and 
is  limited  to  a one-dimensional  treatment.  The  selection  of 
the  Eulerlan  representation  permits  a simple  identification  with 
propellant  locations  along  the  crack  and  allows  for  the  mixing 
of  freshly  burned  reaction  gas  with  gas  from  other  gas  movements. 

It  is  also  convienlent  for  the  large  movements  experienced  by 
some  of  the  gas.  The  one-dimehsional  treatment  is  appropriate 
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because  of  the  extreme  length  of  the  crack  when  compared  to  its 
thickness.  It  is  recognized  that  flow  gradients  will  exist 
across  the  thickness  of  the  crack,  thus  the  value  of  the  variables 
used  in  the  one-dimensional  treatment  are  average  cross-sectional 
values.  The  width  of  the  crack  is  assumed  to  be  much  wider  than 
the  thickness  and  relatively  uniform.  Thus  the  most  significant 
flow  gradients  and  movements  will  occur  in  the  direction  along 
the  length  of  the  crack. 

The  current  model  treats  a crack  of  some  specified  thickness 
variation  at  a reference  time  (defined  as  zero  time) . The  gas  in 
this  crack  at  this  time  is  of  the  same  composition  as  the  reaction 
products  and  exists  in  a rest  state  (no  flow)  at  some  low  pressure 
and  temperature  (say  the  ambient  temperature  of  the  propellant) , 

The  crack  is  connected,  at  the  upstream  end,  to  the  primary  rocket 
motor  chamber  which  is  filled  with  gaseous  reaction  products  which 
are  at  a high  temperature  and  a high  pressure.  The  chamber  state 
is  treated  as  a constant  state  with  respect  to  time  because  of 
the  relatively  short  duration  of  the  events  occurring  within  the 
crack.  This  simple  treatment  can  be  readily  modified  to  accomo- 
date some  other  specified  chamber  state  variation.  The  crack  is 
connected  Instantaneously  to  the  chamber  and  thus  a rapid  Inflow 
of  hot  gas  occurs  in  the  crack.  This  is  the  start  of  thsi  process 
which  may  lead  to  the  unwanted  detonation  of  the  high  energy  pro- 
pellant. It  would  represent  a situation  in  which  an  internal 
crack  propagates  until  it  Intercepts  the  primary  chamber.  A num- 
ber of  other  starting  conditions  could  be  formulated,  however  it 
is  believed  that  tha  above  stimuli  is  somewaht  representative  and 
should  lead  to  conditions  of  interest  in  Identifying  the  mechanisms 
and  interactions  which  may  lead  to  rapid  propellant  burning  within 
the  crack. 

The  downstream  and  of  the  crack  can  be  represented  by  either 
a closed  end  or  a connection  to  a low  pressure  chamber.  In  the 
latter  case  only  gas  outflow  is  permitted  and  tha  pressure  in  the 
low  pressure  chamber  is  held  constant  at  the  level  of  the  initial 
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pressure  within  the  crack.  This  situation  represents  a possible 
connection  to  a relatively  large  debonded  region  which  is  momen- 
tarily isolated  from  the  primary  chamber.  In  this  case  the  crack 
length  is  constant.  If  the  closed  end  crack  condition  is  selected 
it  is  possible  to  permit  the  length  of  the  crack  to  grow  in  some 
arbitrarily  selected  manner.  This  crack  growth  feature  is  not 
based,  as  it  ultimately  should  be,  upon  some  physically  plausible 
crack  propagation  mechanism.  The  incorporation  of  such  a criteria 
together  with  the  possibility  of  the  branching  of  a single  crack 
to  form  a complex  network  of  cracks  represents  some  of  the  ulti- 
mate goals  for  the  improvement  of  the  overall  physical  model. 

B. 2 BASIC  EQUATIONS 

The  equations  of  fluid  dynamics  are  written  with  respect  to 
the  Eulerlan  frame  of  reference  in  which  the  independent  variables 
are  the  distance,  x,  along  the  crack  and  the  time,  t,  These  eq- 
uations represent  the  conservation  of  mass,  momentum,  and  energy. 
The  present  gas  dynamic  model  includes  the  effects  of  inertia, 
wall  friction,  mechanical  wall  response,  and  mass  and  energy 
addition.  Heat  transfer  effects  are  incorporated  into  the 
energy  addition  term.  Furthermore  the  incorporation  of  mechani- 
cal wall  response  effects  appear  at  this  stage  as  wall  movement 
and  flow  area  variations  along  the  crack.  The  mechanical  res- 
ponse of  the  propellant  to  the  imposed  pressure  enviroment  is 
treated  in  Appendix  C. 

The  conservation  laws  can  be  developed  by  considering  a 
control  surface  enclosing  an  element  of  volume  of  length,  dx. 

This  control  surface  is  illustrated  in  Figure  B.l.  The  conserv- 
ation laws  state  that  the  time  rate  of  change  of  the  entity  con- 
sidered, within  the  control  volume  must  equal  the  net  flux  of 
this  entity  through  the  control  surface,  plus  any  related  bound- 
ary contributions  such  as  mass  or  energy  addition,  work  done  at 
the  boundary,  or  boundary  forces.  The  area  of  the  flow  channel, 

A,  is  one  of  dependent  variables  to  be  considered.  However,  in 
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view  of  the  fact  that  the  crack  is  very  wide  compared  to  its 
thickness  or  height, the  width  of  the  flow  channel  can  be  set  equal 
to  unity  and  a corresponding  variable,  the  crack  height,  h,  can 
be  used.  Thus 


A - 1-h  (B-l) 

It  should  be  noted  that  although  this  variable  has  the  dimensions 
of  length,  a dimensional  check  of  the  following  equations  will 
imply  that  it  should  have  the  dimensions  of  length  squared.  This 
is  the  case  when  it  represents  the  area  of  the  flow  channel. 

Tho  three  conservation  equations  are: 

Mass 

+ <»-» 


where  p - gas  density 
u - gas  velocity 

m - mass  addition  per  unit  length  per  unit  time 
The  wall  motion  (velocity)  is  given  by  (3h/3t)  (see  Figure  B-l) 


Momentum 


2H+uaH+Ii£._2rum 

3t  3x  p 3x  ph  pn 


(B-3) 


where  t - wall  shearing  stress.  The  wall  shearing  stress  can  be 
expressed  in  terms  of  the  conventional  pipe  friction  coefficient, 
f,  defined  by 


f - —^7  (B-4) 

%pu 

Furthermore  this  stre3S  acts  in  a direction  opposite  to  the  fluid 
motion.  The  momentum  equation  can  then  be  reformulated  as: 


The  coefficient  f is  assumed  to  be  constant  at  a nominal  value 
because  its  value  and  dependence  are  uncertain  for  this  complex 
flow  environment. 


B5 


Energy. 


9 9 Gm  -1  * 

TF  + U~3ST  + 


9h  p 9h  eTm 


Q pu  3h  p ...  _ 

p 3xvt'“/  pn  " pK  3x  ‘ pH  3t  ph 

where  Q ■ energy  addition  per  unit  length  per  unit  time 
e,j,  ■ specific  total  energy 

The  specific  total  energy  is  given  by 

2 

eT  - e + %u 


(B-6) 


(B-7) 


where  e - specific  internal  energy. 

Equations  (B-2) , (B-5)  and  (B-6)  have  been  arranged  bo  that 
the  terms  which  define  the  contributions  of  the  Bpecial  effects 
are  grouped  on  the  right  hand  aide  of  the  equal  sign,  When  the 
sum  of  the  terms  of  the  right  hand  side  are  set  equal  to  zero 
the  conventional  gas  dynamic  equations  for  the  nonsteady  flow  in 
a constant  area  channel  are  obtained.  These  reduced  equations 
will  be  used  to  establish  the  numerical  method. 


B.3  EQUATION  OF  STATE 

The  equation  of  state  which  has  been  selected  initially  for 
the  propellant  reaction  products  is  that  of  a perfect  gas.  This 
equation  of  state  is  considered  to  be  adequate  for  the  initial 
phases  of  the  investigation,  but  eventually  a better  formulation 
must  be  established. 

The  perfect  gas  representation  is  the  following: 

p - pRT  (B-8) 

where  R - gas  constant 

T » absolute  temperature  of  the  gas. 

Furthermore,  the  specific  heats  at  constant  pressure  and  at  con- 
stant volume  are  both  constant.  The  internal  energy  is  given  as 


• ■ v&r 

where  y - ratio  of  specific  heats. 
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(B-9) 


(B-1Q) 


The  sound  velocity  of  the  gas,  c,  is  given  as 

c -Vie 

p 

The  values  of  the  gas  constant,  the  two  specific  heats  and  the 
ratio  of  specific  heats  are  subject  to  some  uncertainty.  The 
following  values  have  been  selected  for  the  gas  constant  and 
the  ratio  of  specific  heats  for  the  reaction  products 

R - 3228.  (cm°C) 

Y ■ 1 • 2 

B.4  FLOW  CONFIGURATIONS  AND  BOUNDARY  CONDITIONS 

The  solutions  of  the  above  field  equations  are  subject  to 
the  initial  conditions  within  the  crack  and  to  boundary  conditions 
at  both  the  upstream  and  downstream  end  of  the  flow  channel.  The 
initial  conditions  are  those  of  a gas  at  rest  (u(x)  = 0)  and  at 
some  pressure  pQ,  and  temperature,  TQ.  The  initial  height,  hQ(x> 
of  the  channel  is  also  specified,  however  no  initial  wall  motion 
is  permitted  (i.e.  (3ho/3t)=0).  The  propellant  mass  in  the  vicin- 
ity of  the  crack  is  thus  in  mechanical  equilibrium  with  the  ini- 
tial pressure  field  within  the  crack. 

Two  basic  crack  configurations  are  treated.  These  are  illus- 
trated in  Figure  B-2,  They  consist  of  a high  pressure  cavity  con- 
nected to  the  upstream  end  (x-0)  of  a crack  of  length,  Lc,  together 
with  one  of  two  downstream  end  conditions.  One  configuration  has 
a simple  closed  end  while  the  other  configuration  consists  of  a 
connection  to  a low  pressure  cavity.  The  high  pressure  cavity 
conatins  a gas  at  a pressure  pc  and  temperature  Tc,  both  of  which 
are  held  constant  with  respect  to  time.  The  gas  pressure  within 
the  low  pressure  cavity  is  held  constant  at  the  level  of  the  initial 
gas  pressure  within  the  crack  (i.e.  pQ)  and  only  outflow  is  permitted. 

The  boundary  conditions  at  the  high  pressure  cavity  end  will 
depend  upon  whether  the  gas  flow  is  into  (inflow)  of  out  of  the 
crack  (outflow) . These  boundary  conditions  are  illustrated  in 
the  Hodograph  plane  of  Figure  B*>3.  The  cavity  State,  Sc>  ( a rest 
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Figure  B-2  Basic  Crack  Configurations 


state)  Is  the  appropriate  reference  state.  When  the  pressure  is 
low  relative  to  the  cavity  pressure  at  the  downstream  end  of  the 
crack  inflow  may  occur.  If  inflow  does  occur,  the  cavity  gas  will 
flow  into  the  crack  after  first  expanding  isentropically  (at  con- 
stant energy)  in  the  vicinity  of  the  inlet.  This  expansion  process 
will  accelerate  the  gas  and  lower  both  the  pressure  and  the  sound 
velocity.  The  flow  Mach  number,  M,  (M=u/c)  will  increase  until  it 
reaches  a state  which  is  in  equilibrium  with  the  internal  flow 
at  the  boundary  provided  it  does  not  exceed  a value  of  unity.  At 
that  critical  point  no  further  expansion  can  occur  since,  in  effect, 
no  further  information  regarding  any  additional  expansion  can  be 
communicated  back  to  the  cavity.  Any  additional  expansion,  if 
it  can  occur,  will  occur  as  a non-steady  expansion  within  the 
crack.  The  field  equation  will  provide  for  the  solution  of  this 
additional  expansion.  Thus  the  permissible  boundary  conditions  at 
the  downstream  end  of  the  crack  under  inflow  conditions  will  be 
those  states  associated  with  the  energy  line  bounded  by  the  cavity 
State,  Sc,  at  one  end  and  the  sonic  inflow  State  Sfl;L  at  the  other 
end.  The  energy  line  is  defined  by 

CC2  ■ C2  + (^^ip)  U2  (B-ll) 

If  subsonic  outflow  occurs  the  pressure  at  the  upstream  end 
of  the  crack  will  be  equal  to  the  cavity  pressure,  pc  (see  Figure 
B-3) . However,  if  the  flow  becomes  sonic  or  supersonic  then  the 
pressure  can  change  to  any  value  such  that  the  boundary  state 
within  the  crack  lies  in  the  supersonic  outflow  region  bounded 
by  the  sonic  outflow  line.  This  region  is  illustrated  in  Figure  B-3. 

The  same  boundary  conditions  apply  at  the  low  pressure  cavity 
end  of  the  crack  when  this  configuration  is  used,  however  the 
present  model  does  not  permit  inflow  to  occur.  Some  model  modi- 
fications in  this  area  may  ultimately  be  made  however,  they  will 
be  influenced  by  the  physical  model  associated  with  this  low 
pressure  cavity,  such  as  for  example  a cavity  filling  process  with 
a subsequent  buildup  of  pressure  and  temperature  within  the  finite 
volume  cavity. 
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The  boundary  condition  for  the  closed  end  is  a simple  one  of 
no  flow  (i.e.  u(L  )=0).  The  current  model  does  permit  the  crack 

C 

length  to  be  extended  arbitrarily  whenever  this  configuration  is 
used.  Whenever  the  crack  is  extended, the  new  portion  is  filled 
with  a rest  gas  at  the  initial  pressure  and  temperature  (p0»TQ) 
and  the  closed  end  boundary  condition  is  applied  to  the  new  closed 
end  location.  This  crack  extension  is,  in  effect,  an  instantaneous 
crack  extension. 

B.5  NUMERICAL  PROCEDURE 

The  numerical  solution  technique  employed  in  the  gas  dynamic 
model  of  the  present  study  is,  in  its  reduced  form,  a rather  con- 
ventional Eulerian  method  of  the  FLIC  (fluid  in  cell)  type.  This 
method  was  originally  developed  at  the  Los  Alamos  Scientific 
Laboratory. 

The  numerical  solution  of  the  foregoing  equations  proceeds 
from  the  initial  conditions  specified  in  a forward  stepping  time 
wise  manner  subject  to  the  numerical  solution  of  the  field  equa- 
tions, auxilliary  timewise  inputs,  and  the  appropriate  boundary 
conditions.  An  artificial  viscous  pressure  term,  q,  is  added  to 
the  thermodynamic  pressure  during  the  computations  when  the  flow  is 
subsonic  and  the  compression  rate  is  positive.  This  contributes 
to  the  suppression  of  flow  discontinuities  and  to  the  computational 
stability  in  regions  of  subsonic  flow. 

B.5.1  The  Computing  Mesh 

A one-dimensional  mesh  of  uniform  length  Ax,  cells  is  estab- 
lished. Each  cell  is  Identified  by  an  lndice,  i,  corresponding  to 
its  center.  Thus  the  boundaries  of  the  ith  cell  are  at  the  loca- 
tion i+%.  Increasing  1 corresponds  to  increasing  x.  The  selection 
of  uniform  cell  length  is  an  initial  convienience . If  subsequent 
results  indicate  that  one  portion  of  the  crack  requires  substan- 
tially greater  resolution  in  the  space  variable,  x,  than  do  other 
portions  then  variable  length  cells  can  be  introduced  with  little 
additional  complications. 

BU 


The  values  of  the  dependent  varlahles  are  associated  with 
the  cell  centers  with  the  exception  of  the  pseudoviscous  terms 
which  are  computed  at  the  cell  boundaries.  Assuming  that  all 
properties  are  known  for  each  cell  at  some  time  tn,  the  computa- 
tional procedure  is  to  determine  the  state  In  each  cell  at  a 
later  time  tn+*  - tn  + At.  The  time  step  At  is  restricted  in 
magnitude  by  conditions  required  for  stability  of  the  computations. 
The  resulting  space  time  grid  network  is  illustrated  in  Figure  B-4. 


B . 5 . 2 Stability  Formulation 

There  are  several  stability  type  conditions  that  other  inves- 
tigators have  found  applicable  with  this  type  numerical  technique. 
One  such  restriction  is  that  (u  avAt) /Ax  < 1 . If  fluid  particles 
were  explicitly  treated  in  the  computation,  this  criterion  would 
prevent  a particle  from  crossing  a complete  cell  during  one  time 
step.  In  the  pure  Eulerian  scheme  the  Interpretation  is  made 
that  the  transport  terms  are  calculated  more  accurately  with  cor- 
responding improved  averaging  of  the  numerical  fluctuations.  A 
second  restriction  is  the  familiar  Courant  condition,  cAt/Ax<l, 
which  limits  the  propagation  of  fluctuations  in  the  subsonic 
flow  regions.  Another  source  of  instability  is  the  computation  of 
negative  internal  energies  in  the  cell  due  primarily  to  the  treat- 
ment of  boundary  conditions,  but  also  from  inherent  numerical  fluc- 
tuations. For  a fluid  obeying  the  ideal  gas  relation,  the  criterion 
that  prevents  negative  internal  energies  is  that  4(y-l)  umaxAt/,Ax<^' 
These  stability  criteria  can  be  written  as 


At  < Min 
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(B-13) 

(B-14) 


In  the  last  equation  a value  of  1.4  for  y was  used.  A simple 
stability  criteria  which  satisfies  all  of  th«  foregoing  require- 
ments end  which  is  used  in  the  calculations  is 
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where  the  factor  % is  chosen  to  reflect  the  typical  severity  with 
which  the  inequalities  (B-12)  to  (B-14)  are  generally  applied. 

The  introduction  of  special  effects  may  have  an  impact  upon  the 
stability  of  the  calculations,  therefore  additional  time  step 
criteria  may  have  to  be  introduced. 

B . 5 . 3 Artificial  Viscosity  Term 

Artificial  viscous  pressure  terms  are  computed  at  time  tn 
at  those  cell  boundaries  acrocs  which  a negative  velocity  compo- 
nent gradient  exists,  provided  that  the  local  Mach  number  does 
not  exceed  a prescribed  value  defined  by  a constant,  A . The 
terms  are  proportional  to  the  velocity  component  differences 
between  adjoining  cells  at  the  boundary  where  it  is  applied. 
Specifically, 
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provided  - u1+1  > 0 
and  Aqu2i+3k  - C2^  < 0 

If  either  of  the  indicated  flow  conditions  is  not  satisfied,  the 
viscosity  term  is  set  equal  to  zero.  The  values  of  Ag  and  B used 
in  the  calculations  were  100  and  0.5  respectively.  The  subsonic 
limit  correspond  to  a Mach  number  of  0.1. 

B.5.4  Intermediate  Values  of  the  Flow  Variables 

At  this  stage  of  the  development  of  the  numerical  method  only 
the  reduced  form  of  the  field  equations  will  be  created.  These 
equations  are  rewritten  in  the  following  form. 


ft  + ^ <pu)  * 0 
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where  intermediate  values  of  the  velocity  and  the  internal  energy 
(denoted  by  the  circumflex  'v)  are  computed  from  their  local  rates 
of  change  neglecting  the  transport  terms  in  the  field  equations. 
The  gradient  and  divergence  terms  in  the  momentum  and  energy  equa 
tions  are  calculated  from  numerical  approximations  to  the  equiv- 
alent surface  integrals  of  pressure  and  velocity  over  the  cell 
boundaries.  The  resulting  difference  equations  are 

“i  - ui  - ^ [ wu  - <b-2°> 

and 


where  the  macrons  used  in  Equation  (B-21)  represent  the  average 
(un  + un)/2  at  the  indicated  subscript.  These  forms  provide 
significant  advantages  in  the  development  of  numerical  solution 
algorithms . 

An  Intermediate  value  of  the  total  energy  eT  is  computed 
from  its  definition,  Equation  (B-7) , using  the  Intermediate 
values  of  internal  energy  and  velocity  Just  calculated. 

B.5.5  New  Values  of  the  Flow  Variables 

Final  values  of  the  flow  variables  at  the  new  time  tn+* 
are  calculated  next  for  the  call  centers  by  including  the  trans- 
port effects.  The  procedure,  in  principle,  is  to  associate  with 
the  intermediate  flow  field  a transport  of  mass,  momentum  and 
total  energy,  defined  in  terms  of  their  intermediate  values, 
across  the  cell  boundaries.  The  technique  of  donor  cell  dif- 
ferencing is  used  in  which  the  transport  of  these  quantities  is 
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calculated  from  the  donor  cell,  t.e.,  the  cell  from  which  the 
fluid  is  flowing.  This  technique  prevents  cells > from  developing 
negative  densities  due  to  numerical  fluctuations  and  provides 
good  stability  properties  in  the  far  subsonic  regions  of  the  flow. 


The  mass  flow  terms  are  calculated  from  the  numerical  approx- 
imation of  the  surface  integral  equivalent  to  the  divergent  term 
appearing  in  Equation  (B-2) , resulting  in 

“J+i  ■ (p1H)D^+%  4t  (°-22) 

The  superscript  D is  used  to  indicate  that  the  density  is  assigned 
the  value  according  to  the  donor  cell  requirement,  which  may  be 
written  as : 
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/ n vD 

(pi+%> 


K+k  > 0 


n j-n 
kpi+l*  u 


i+% 


< 0 


It  is  now  appropriate  to  rewrite 
equations,  Equations  (B-17)  to  (B-19) 
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the  momentum  and  energy 
as 


and 
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where  the  reduced  nonconvective  equation,  (Equations  (B-17)  to 
(B-19))  that  are  used  to  define  the  intermediate  flow  variables  that 
have  been  used,  The  transport  terms  appearing  in  Equations  (B-24) 
and  (B-25)  are  evaluated  from  numerical  approximate  to  their 
equivalent  surface  integral  representations.  Next,  applying 
first  order  finite-differences  tosll  temporal  derivatives,  the 
following  general  difference  formulation  of  the  reduced  field 
Equations  (B-17) , (B-24)  and  (B-25)  is  obtained. 
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where  the  variable  F is  used  to  represent  in  turn  the  various 
flow  variables  according  to 


F 


i 


(B-27) 


The  identity  value  of  F can  be  interpreted  as  a mass  flow  which 
transports  only  itself,  thereby  conserving  mass  and  allowing 
Equation  (B-26)  to  be  used  in  the  calculations  of  the  new  cell 
density  p™  . The  donor  cell  requirement  on  F is  given  by 
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Substituting  the  mass  flow  equation  (B-22)  into  Equation 
(B-26)  yields  the  following  calculational  form  for  which  extremely 
efficient  numerical  solution  algorithms  may  be  developed: 
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Once  Equation  (B-29)  has  been  evaluated  for  the  density, 
velocity  and  total  energy,  the  new  values  of  internal  energy, 
pressure,  sound  velocity  and  temperature  are  readily  calculated 
from  Equations  (B-7) , (B-9) , (B-10) , and  (B-8) , respectively. 


B.5.6  Treatment  of  the  Special  Terms 

The  previous  description  of  the  numerical  procedure  dealt 
with  the  reduced  form  of  the  basic  equations  (Equation  (B-2) , 

(B-5)  and  (B-6)).  The  special  terms  of  these  equations  can  be 
incorporated  into  the  numerical  procedure  in  a relatively  simple 
manner,  First  consider  the  friction  term  of  the  momentum  equation. 
This  term  contributes  to  the  local  rate  of  change  and  can  therefore 
be  incorporated  into  the  calculation  of  u.  Thus  Equation  (B-19) 
can  be  expanded  as 


B17 


(B-30) 


1M  - - I l£  . fuJul 

at  p ax  n 
and  changes  Equation  (B-2Q)  to 


^n 

u 


- 


At 

pjAx 


[(p+q)J+%  - (p+q>i.%]-  At 


fujlujl 

“ t_n 


(B-31) 


The  mass  and  energy  equations  can  be  rewritten  in  the  form. 

ft  + K St  <*“h>  - f * f <B-32> 


where  W *■  ah/ at  » wall  velocity  and 
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The  intermediate  value  of  the  internal  energy  is  now  computed  from: 
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Therefore  Equation  (B-21)  becomes 
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The  final  values  of  the  flow  variable  at  the  new  time  t " " are 
then  calculated  from  considerations  of  the  following  expanded 
equations  for  momentum  and  energy 
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and  Equation  (B-32) 
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The  general  solution  equation,  Equation  (B-29)  is  no  longer  appli- 
cable because  of  the  various  special  terms.  Rather  three,  some- 
what similar  relations  are  needed.  These  are  the  following. 
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where  is  defined  by  Equations  (B-27)  and  (B-28) . 

and 


n+1/ 

Pi 


.n+1  _ 

;i 


P,  <e 


T^i 


At 

Ax 


k 


Fi-%  Ui-% 


- p2+%  Kh,  V*}  + [«i  - <ST>i“i]  <B-40> 


B.5.7  Boundary  Conditions 

The  boundary  conditions  imposed  upon  the  two  ends  of  the 
flow  channel  are  accommodated  by  the  use  of  one  dummy  cell  at 
each  end.  The  first  upstream  cell  of  the  flow  channel  is  desig- 
nated by  i-2  and  the  associated  dummy  cell  is  designated  by  i*»l. 

If  in  cells  are  used  to  represent  the  flow  channel  then  the  last 
real  cell  is  designated  by  i^-in+l.  The  downstream  dummy  cell 
is  designated  i^i^+1.  Thus  a total  of  ic  cells  are  used. 

The  cell  i-1  represents  the  boundary  cell  for  the  high  pres- 
sure cavity  and  some  values  of  the  variable  in  this  cell  will 
depend  upon  whether  inflow  or  outflow  conditions  exist.  This 
condition  is  determined  by  examining  the  value  of  U2  at  any  time 
tn.  For  outflow  « pc  provided  the  flow  is  subsonic,  otherwise 
Pj_  ■ p? . This  latter  relation  provides  better  numerical  uncoup- 
ling from  the  cavity  under  sonic  and  supersonic  outflow  conditions. 
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For  inflow  p-^  is  determined  from  the  energy  equation  in  the 
following  form 


Pi  = Pc 
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where  - U2  provided  that  its  value  does  not  exceed  the  sonic 
value,  u . This  critical  value  is 


u. 


m 


c 

c 


(B-42) 


The  cavity  sound  velocity,  cc,  is  determined  from  Equations 
(B-8)  and  (B-10) . The  velocities  u and  u in  cell  1 for  both  in- 
flow and  outflow  are  set  equal  to  their  corresponding  values  of 
cell  2,  via. 


(B-43) 


provided  that  they  do  not  exceed  the  algebraic  value  of  um> 

This  limiting  requirement  can  only  be  experienced  during  inflow, 
The  density  and  total  energy  (e^)^  are  required  for  inflow. 
They  are  given  by  the  following  relations. 
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The  boundary  conditions  are  the  downstream  end  are  dependent 
upon  the  type  of  end  condition  being  used.  For  the  closed  end 
configuration  the  particle  velocity  at  the  end  of  the  channel 
must  vanish.  This  is  accomplished  by  setting  the  intermediate 
values  of  the  velocities  ft  and  u in  cell  ic  equal  to  the  negative 
of  the  values  of  their  corresponding  values  in  cell  1^,  viz. 
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Uic  " -uib 
uic  " "uib 

In  addition  the  required  state  variables  of  cell  ic  are  set  equal 
to  the  corresponding  values  of  cell  ife  to,  in  effect,  create  a 
plane  of  symmetry,  specifically 

Pic  " pib 

Ptc  - Pib  (B-47) 

^eT^ic  " ^eT^ib 

The  downstream  boundary  conditions  for  the  low  pressure  cavity 
are  considerably  simpler  than  those  for  the  high  pressure  cavity 
because  only  outflow  is  permitted  to  occur.  The  following  values 
are  needed. 

pic  " p0  if  Mib  < °'  otherwise  pic  - PjLb 
“ic"Ab  <*-48) 

“ *ib 

Where  pQ  is  the  initial  pressure  of  the  gas  within  the  crack. 

B.5.8  Calculation  Method 

The  basic  calculation  procedure  is  to  evaluate  Equations 
(B-16) , (B-31)  , (B-35) , (B-38),  (B-39) , (B-40) , (B-9) , (B-7) , 

(B-10) , and  (B-8)  and  the  appropriate  auxiliary  equations  for  each 
cell  of  the  computing  mesh  (cells  i-2  to  ib)  during  each  time  step. 
The  numerical  results  computed  at  the  end  of  any  time  step,  together 
with  a new  set  of  boundary  conditions,  then  serve  as  the  initial 
conditions  for  the  next  time  step,  and  in  this  manner  the  numerical 
solution  ia  advanced  in  time. 
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a. 6 CODE  VERIFICATION 


A computer  program  was  developed  to  implement  the  preceding 
gas  dynamic  model.  This  computer  program  exists  (1)  in  the  form 
of  a main  program  with  which  to  examine  some  aspects  of  the  gas 
dynamic  phenomena  and  to  serve  as  a vehicle  for  further  develop- 
ment, and  (2)  in  the  form  of  a subroutine  of  the  overall  propel- 
lant detonation  computer  model.  These  codes  were  programmed  in 
FORTRAN  IV  user-oriented  source  language  and  have  been  success- 
fully executed  on  a UNIVAC  1108,  a large  scientific  computer. 

The  gas  dynamic  code  has  been  exercised  to  evaluate  its  ade- 
quacy for  the  subject  problem  area.  Current  results  indicate  that 
reasonable  resolution  can  be  obtained  for  a class  of  simple  flow 
problems  if  approximately  20  or  more  cells  are  used  to  define  the 
crack  length.  Comparisons  have  been  made  with  a number  of  analy- 
tic solutions  for  hoth  transient  and  steady  state  (late  time) 
flows  and  indicate  that  the  present  form  is  adequate  except  for 
the  treatment  of  the  variable  area  and  thus  moving  wall  features. 
It  is  not  clear  at  the  present  time  whether  the  difficulty  is  due 
to  the  limited  number  of  cells  used  to  define  a large  change  in 
area,  a time  step  requirement,  or  to  some  other  cause.  The  source 
of  the  error  must  be  examined  in  more  detail  and  the  gas  dynamic 
code  will  then  be  modified  accordingly.  This  problem  is  not 
considered  to  be  a limiting  problem. 
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APPENDIX  C 

MECHANICAL  RESPONSE  MODELS 


C.l  INTRODUCTION 

The  flow  of  high  pressure  gases  into  cracks  within  and  adja- 
cent to  the  propellant  mass  and  any  subsequent  burning  of  these 
propellant  surfaces  will  create  a time  varying  pressure  environ- 
ment with  the  cracks.  The  surrounding  propellant  mass,  which  was 
originally  in  mechanical  equilibrium  with  a low  pressure  distribu- 
tion within  the  crack  will  respond  to  these  new  mechanical  loads. 

It  is  the  purpose  of  this  appendix  to  describe  several  mechanical 
response  models  which  have  been  used  in  conjunction  with  the  gas 
dynamic  model  described  in  Appendix  B.  One  response  model  treats 
the  early  phases  of  the  mechanical  response  when  stress  waves 
radiate  from  the  crack  and  the  crack  is  in  a growth  phase.  The 
second  model  tests,  in  a simple  fashion,  the  Interaction  of  these 
radiating  waves  with  the  confining  shell  structure  and  their  sub- 
sequent interaction  with  the  crack.  Under  these  conditions  the 
crack  growth  will  generally  be  arrested  and  may  ultimately  lead 
to  partial  or  complete  crack  collapse. 

C.l  MECHANICAL  PROPERTIES  OF  THE  PROPELLANT 

The  propellant  response  to  the  transient  pressure  loads  within 
a crack  will  be  that  associated  with  the  generation  and  propaga- 
tions of  stress  waves  in  the  propellant  mass.  The  mechanical 
properties  of  most  propellant  materials  are  quite  complex  in  that 
they  are  generally  nonlinear,  rate  sensitive,  and  hysteretic  in 
nature.  Any  simple  mechanical  response  model  cannot  deal  affec- 
tively with  these  types  of  complex  behavior  characteristics; however 
their  Impact  upon  the  subject  problem  is  not  considered  to  be 
Important.  Rather  it  should  be  sufficient  to  model  in  some  simple 
way  the  gross  compressibility  characteristics  of  the  material. 

Thus  the  material  will  be  viewed  Initially  as  a simple  linear  iso- 
tropic elastic  material  which  does  not  yield  or  fall  in  any  manner 
under  the  imposed  stresses. 
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An  elastic  material  will  support  two  basic  types  of  stress 
waves,  i.e.,  dilatation  and  shear  waves.  These  waves  will  pro- 
pagate at  the  following  speeds. 

Dilatation  Wave  Speed,  cj 

cd  "Vl  Vu-VyKio  <c“l> 

Shear  Wave  Speed,  cs 

cs  "VTVztctt 7 CC‘2) 

where  E - Young's  modulus 
t - Poisson' 8 ratio 
p - density 

The  ratio  of  these  two  wave  speeds  is  dependent  only  upon  Poisson's 
ratio,  viss 


1-2t 


(C-3) 


The  influence  of  Poisson's  ratio  is  not  great  until  it  exceeds  a 
value  of  about  .45  and  then  it  primarily  affects  the  dilitation 
wave  speed.  The  effective  wave  speed  c*  • ✓’E/p  is  the  most  signi- 
ficant parameter.  Since  Poisson's  ratio  may  be  in  the  broad  range 
of  from  0.2  to  0.4  the  dilitation  wave  speed  will  be  approximately 
1.2  C*  while  and  the  shear  wave  speed  will  be  approximately  0.6  c*. 
Thus  the  wave  speed  ratio  will  be  approximately  0.5.  The  density 
of  the  propellants  are  well  known  and  a nominal  value  of  1.5  g/cm^ 
is  used  in  the  mechanical  response  models.  A nominal  value  of  250 
cm/ms  was  selected  for  the  effective  wave  speed.  This  corresponds 
to  a value  of  approximately  IQ®  pal  for  Young's  modulus  and  yields 
dilitation  and  shear  wave  speeds  of  300  cm/ms  (approximately  10,000 
fps)  and  150  cm/ms  respectively. 
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C, 3 INITIAL  RESPONSE  MODEL 


During  the  entry  phase  of  the  gas  flow  into  the  crack  a shock 
wave  is  generated  which  propagates  at  a speed  (depending  upon  a 
number  of  parameters)  of  approximately  150  cm/ms.  The  pressure 
distribution  behind  this  shock  wave  is  relatively  uniform.  Weak 
disturbances  (i.e.,  soundwaves)  propagate  at  the  local  sound 
speed  relative  to  the  gas.  The  sound  speed  is  generally  in  the 
approximate  range  of  from  50  to  100  cm/ms  while  the  particle  veloc- 
ity is  in  the  range  of  from  -50  to  100  cm/ms.  Thus  these  disturb- 
ances will  propagate,  on  the  average,  at  an  absolute  velocity  of 
about  100  cm/ms.  These  considerations  lead  us  to  the  following 
approximate  inequality , 

| (y+c)  | < U dcg<c(j  (C-4) 

where  U - shock  velocity  in  the  gas 

y ■ particle  velocity  of  the  gas 
c «*  sound  velocity  of  the  gas 

Thus  the  wave  system  illustrated  in  Figure  8 of  the  main  text 
should  exist  initially.  While  this  wave  system  is  quite  complex 
many  of  the  waves  should  be  weak  and  need  not  be  considered.  This 
is  especially  true  for  the  outrunning  waves . The  primary  motion 
of  the  propellant  mass  will  be  in  a direction  normal  to  the  crack, 
causing  the  crack  height  to  increase  locally.  For  this  reason, 
a very  simple,  one  dimension  wave  propagation  (plane  strain)  model 
has  been  selected  with  which  to  define  the  response  of  the  propel- 
lant adjacent  to  the  crack.  The  momentum  equation  for  this  case 
is  of  the  form 

At  - (pcd)Au  (C-5) 

where  Ax  - change  in  stress  at  the  crack  boundary 

Ay  • change  in  velocity  at  the  crack  boundary 

The  change  in  the  stress  at  a given  location  along  the  crack  is 
identical  to  the  change  in  the  corresponding  gas  pressure.  Thus 
the  wall  velocity  W used  in  the  gas  dynamic  model  is  given  for 
the  ith  cell  as 


«1  - (Pcd>(PrP0) 
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where  pc^  - shock  impedance  (.*  450  bars-ms/cm) 

- gas  pressure 
Po  ■ initial  gas  pressure 

This  expression  was  modified  to  include  some  contribution  from 
the  adjacent  cells  (i-1,  and  i+l)  . In  this  manner  some  influence 
of  the  spreading  nature  of  the  propagation  within  the  propellant 
mass  was  included,  at  least  in  a crude  manner. 

(3.4  INFLUENCE  OF  REMOTE  BOUNDARIES 

The  stress  wave  which  radiates  away  from  the  crack  will 
eventially  react  with  the  remote  boundaries  of  the  propellant 
mass,  reflect  and  return,  in  some  modified  form,  to  the  crack. 

These  reflected  waves  will  also  influence  the  response  crack  walls. 

This  type  of  wave  system  is  illustrated  in  Figure  9(a)  of  the 
main  text.  The  nature  of  reflected  wave  system  will  be  very 
complex  and  varied.  For  this  reason  a number  of  idealized  re- 
flection models  will  have  to  be  established  to  define  the  con- 
figurations and  conditions  of  interest. 

As  an  initial  attempt  to  treat  the  influence  of  the  remote 
boundaries  of  the  propellant  mass  one  idealized  model  was  estab- 
lished. This  model  is  illustrated  in  Figure  9(b).  It  is  designed 
to  introduce  the  ultimate  confining  influence  of  a cased  mass  of 
propellant.  A radiative  stress  field,  Apr,  which  is  the  integra- 
ted or  averaged  value  of  the  current  local  stress  field  along 
the  entire  crack  is  defined  and  assumed  to  propagate,  under  the 
conditions  of  plane  strain,  into  the  propellant  mass.  These  stress 
waves  interact  with  the  remote  boundary  after  a delay  time  corres- 
ponding to  a boundary/crack  separation  distance,  Le>  The  case 
is  treated  as  a spring  supported  mass  characterized  by  its  inertia 
(weight  per  unit  area)  and  stiffness  (breathing  mode  of  the  case). 

The  reflected  wave  system,  Aps  is  evaluated  from  this  case  boundary 
interaction  and  applied  uniformly  along  the  crack  after  the  appropriate 
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delay  period.  This  reflected  wave  system  then  interacts  locally 
with  the  crack  boundary  such  that  the  wall  velocity  changes  ac- 
cording to  the  following  momentum  consideration 

AWt  = -2Apg/(pcs)  (C-7) 
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APPENDIX  D 

TECHNIQUE  FOR  PREDICTING  PROPELLANT  TEMPERATURES 
D.l  INTRODUCTION 

In  this  appendix  we  shall  discuss  the  techniques  used  to 
compute  transient  temperature  rises  within  regressing  propellants 
produced  by  known  time-dependent  thermal  fluxes.  Topics  covered 
in  order  of  presentation  are 

• prediction  of  transient  temperature  rises 

• selection  of  time  steps 

• means  used  to  expedite  computations 

• validation  of  the  technique  by  comparison 
with  an  exact  solution 

Means  for  predicting  the  transient  velocities  of  the  regressing 
boundary  are  discussed  in  Appendix  A. 

In  essence  the  method  involves  a series  of  constant  thermal 
fluxes  applied  at  various  locations  of  a regressing  boundary  with- 
in a semi-infinite  body  during  selected  time  intervals.  These 
fluxes  include  incident  fluxes  along  with  appropriate  combinations 
of  conductive  fluxes  produced  at  other  depths  and  times  by  the 
incident  fluxes.  The  net  result  in  a series  of  independent  heat 
conduction  problems  which  are  then  solved  and  summed  to  arrive 
at  the  temperature  rises  within  the  remaining  propellant. 
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D. 2 PROPELLANT  TEMPERATURES 


Prediction  of  temperatures  beneath  accelerating  boundaries 
have  been  handled  by  a variety  of  techniques  (37-39) . Of  all  the 
methods  examined,  finite-difference  techniques  are  the  most  versa- 
tile  and  appropriate.  Unfortunately,  finite-difference  techniques 
require  considerable  computer  storage  and  execution  times  due  to 
the  large  number  of  spacial  and  temporal  elements  needed  to  accomo- 
date rapid  changes  of  velocities  and  temperatures.  This  is  true 
with  implicit,  explicit  or  implicit-explicit  methods  — all  of 
which  require  very  long  computational  times  to  treat  highly  dynamic 
problems . 

As  a result  a novel  technique  was  developed  by  modifying  a 
closed-form  solution  of  transient  heat  conduction  within  a semi- 
infinite body  with  a stationary  boundary  exposed  to  a constant  flux. 
Here  a fixed  coordinate  system  is  used  within  which  the  boundary 
is  allowed  to  regress  in  a discrete  fashion  during  selected  time 
intervals.  During  each  time  interval,  the  Incident  flux  is  held 
constant  corresponding  to  a step-wise  approximation  of  the  time- 
dependent  flux  q(t)  at  the  propellant  boundary.  By  determining 
the  rates  of  heat  conduction  occurring  at  various  depths  and  times, 
one  arrives  at  a series  of  simple  heat  conduction  problems  from 
which  the  temperatures  rises  within  the  moving  propellant  are 
closely  approximated. 

To  illustrate  the  process  the  reader  is  referred  to  Fig.  D-l. 

In  this  figure  the  boundary  is  considered  to  commence  to  move 
starting  at  tj_^  and  to  arrive  at  Xj  at  the  end  of  the  time  in- 
terval Atj , Incident  fluxes  And  mean  depths  are  represented  by 
qt  and  *1> 

Once  the  boundary  is  displaced,  fluxes  lower  than  q^  are 
used  to  negate  undesireable  conductive  fluxes  arising  in  the  com- 
putational scheme  from  prior  incident  fluxes. 
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The  technique  starts  with  the  heat  conduction  equations  for 
temperature  rises  T(x,t)  produced  within  a semi-infinite  body 
exposed  to  a constant  flux  q^  applied  at  some  depth  x^  for  times 
t > t^_^.  These  equations  are  presented  below. 
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where  a,K  are  the  thermal  diffusivity  and  conductivity  of  the 
propellant. 

Assumption  of  a semi-infinite  body  is  a valid  one  since  heat 
penetration  is  relatively  shallow  compared  to  propellant  depths. 
Solution  of  the  above  problem  with  x^t^ _^«0  is  given  in  Ref.  36. 
By  simple  transformation  of  this  solution  one  arrives  at  the  fol- 
lowing solution  of  the  above  equations. 

X M X 

T(x,t)-C'q<  A-t.  , ierfc  1 (D-5) 

2 /oCt-t^) 

for  x > x ^ , t > t^_^ 

where  G ' -2/ /KpC  and  K,p,C  represent  the  thermal  conductivity,  den- 
sity and  specific  heat  of  the  propellant,  respectively. 

D.2.1  Stationary  boundaries 

For  this  case  x^  is  zero.  In  view  of  the  boundary  condition 
of  Eq.  D-2,  it  is  possible  to  compute  the  temperature  rises  pro- 
duced by  each  q^  at  x^-0  separately,  and  sum  the  resultant  temper- 
atures. To  demonstrate  this  fact  we  shall  replace  Eq.  D-2  with 
the  following  boundary  conditions  associated  with  each  of  a series 
of  fluxes  applied  in  a sequential  manner. 
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Applying  Duhamel's  principal  to  Eq.  D-5  to  accommodate  the  boun- 
dary conditions  of  Eq.  D-6  for  any  single  q^  yields 
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Assuming  the  temperature  rises  produced  by  each  are  addi- 
tive, then  the  temperature  rises  T(x,t)  produced  by  the  first  k 
fluxes  q^,  where  k<  1-1,  is 
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, for  tk_1<t<tk 


To  prove  that  Eq.  D-8  satisfies  the  boundary  conditions  of 
Eq.  D-6,  let  us  multiply  Eq.  D-8  by  -K  and  differentiate  with  re- 
spect to  x to  arrive  at  the  flux  given  below. 
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For  the  first  time  interval  k**l.  And  since  to=0,  Eq.  D-9 
reduces  to 

q,  erfc  — - — , 0<t<t,  (D-10) 
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At  X“0,  it  is  clear  that  the  flux  is  q^.  During  the  second  time 
interval  wherein  k-2,  Eq.  D-9  reduces  to 
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for  t^<t<t2 

At  x*0,  the  flux  is  q2  during  the  above  time  interval. 


By  continuing  the  process  for  the  first  j-1  time  intervals, 
the  boundary  conditions  cited  by  Eq.  D-6  can  be  shown  to  be  satis- 
fied for  each  of  the  fluxes  cj^.  And  since  Eq.  D-8  satisfies  Eqs. 
D-l , D-3,  and  D-4  on  a term  by  term  basis,  it  correctly  describes 
the  temperature  rises  produced  by  the  first  j-1  fluxes  q^. 


0-2.2  Moving  Boundaries 

Once  the  boundary  starts  to  move,  one  must  commence  to  dis- 
count the  most  recent  fluxes  produced  at  the  displaced  boundary 
by  prior  q^.  To  illustrate  the  procedure,  we  shall  first  treat 
the  problem  following  the  initial  displacement  of  the  boundary 
to  x, . Then  we  shall  generalize  the  results  for  multiple 
displacements . 

For  the  Initial  displacement  we  shall  assume  a solution  of 
Eqs.  D-l,  D-3,  D-4  and  D-6  of  identical  form  as  Eq.  D-8,  namely 
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where  x^  represents  the  effective  depth  at  which  the  are  applied 
and  is  given  by 

x^  - + C2*Axi  (D-13) 

Within  the  summation  sign  of  Eq.  D--12,  each  of  the  x^  are  aero  and 
^i“qi  ^see  D"^  • F°r  Pr^-or  times,  Eq.  D-12  reduces  to  Eq.  D-8. 

Here  the  fluxes  q^  are  applied  at  given  depths  within  the 
displacements  x^  to  which  they  are  associated.  This  means  that 
the  constant  C2  has  a positive  value  less  than  1.  The  consequence 
of  using  various  C2  values  is  presented  in  Section  5 of  this  appen- 
dix. In  this  regard  it  was  found  that  a Cg  value  of  about  0.3 
yields  the  most  accurate  approximations  of  the  temperature  profiles 
and  heat  content  within  the  receding  material. 

Equation  D-12  obviously  satisfies  Eqs.  D-l,  D-3  and  D-4.  To 
satisfy  the  boundary  condition  of  Eq . D-6  associated  with  the  most 
recent  time  interval  Atj,  we  shall  first  multiply  Eq.  D-12  by  ~K 
and  then  differentiate  the  result  with  respect  to  x.  Setting 
t-tj  and  x-Xj  yields  the  following  approximation  for  the  flux  qj : 

1-1  _ _ _ _ 

x,-x<  x,-x.  \ 

erfc  — ~ J— - erfc  — ^ ■■  ■ 1+  q.  (D-1A) 

, 2/a(tj-ti“1)  2A(tj-ti)  j ■> 

Here  the  approximate  expression  within  the  summation  sign 
represents  relatively  small  undesirable  conductive  fluxes  arising 
in  the  computational  scheme  from  prior  q^  at  x^  during  At^ . 
Recognizing  that  the  undesirable  fluxes  depend  upon  time,  a more 
exact  correction  can  be  achieved  by  integrating  the  fluxes  over 
the  most  recent  time  interval  tj-t^j-Atj . This  can  be  accom- 
plished by  replacing  each  of  the  1th  terms  within  the  summation 
sign  of  Eq.  D-14  with 
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Letting  8 equal  the  argument  of  either  erfc  function  of 
Eq.  D-15  and  performing  the  integrations  by  parts  yields 
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Replacing  the  terms  within  the  summation  sign  of  Eq.  D-14  by  the 
expression  given  by  Eq.  D-16  and  solving  for  q^  yields 
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To  speed  computer  execution  time  the  function  y(8)  may  be 
approximated  by  the  following  series 
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(D-20) 


Here  sufficient  terms  have  been  taken  in  the  series  expansion 
so  that  the  above  approximation  agrees  to  within  one  percent  of 
that  presented  by  Eq,  D-18.  Errors  Introduced  into  the  computed 
temperature  changes  by  the  above  approximation  for  y(8)  are  less 
than  one  percent. 

Substituting  the  expression  for  q,  given  by  Eq . D-19  into 
Eq.  D-12  yields  the  temperature  rises  following  the  initial  dis- 
placement Ax.  of  the  boundary.  By  continuing  the  procedure  one 
arrives  at  the  method  for  multipla  displacements  given  in  the 
next  section. 


D.2.3  Summary 


In  general,  the  temperature  rises  are  given  by 
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for  x>xn<  The  term  qn  is  given  by  the  following  recursion  formula: 
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D . 3 SELECTION  OF  TIME  INTERVALS 

An  important  benefit  of  this  method  is  its  capability  of  using 
relatively  large  time  intervals  compared  to  finite-difference 
methods.  There  are  no  problems  of  instability.  The  principal 
restriction  on  the  time  intervals  is  that  they  must  be  kept  suffi- 
ciently small  to  accurately  describe  the  incident  thermal  fluxes 
and  prevent  flux  variations  across  the  displacements  Ax^  from 
becoming  excessively  large. 

During  the  initial  heating  period  before  the  boundary  commences 
to  move,  the  time  intervals  are  selected  solely  on  the  basis  of 
yielding  a reasonable  step-wise  approximation  of  the  time-dependent 
changes  in  the  flux.  Once  the  boundary  commences  to  move,  the 
time  intervals  are  computed  as  follows: 

Att  - a(C1q1_1/ CV1-1q±))2  (D-23) 

2 

where  a - thermal  diffusivity  in  cm  /sec 

qt  ■ thermal  flux  in  cal/cm^-sec  during  the  time 
1 interval  At.^  in  sec 

" constant  usually  equal  to  0.3  or  less 
Vt-i  « previous  velocity  in  cm/sec 

This  expression  limits  the  most  recent  arguments  of  the  ierfc  and 
erfc  functions  to  values  of  about  C^/2.  Accuracy  will  increase 
as  the  constant  is  decreased  since  such  will  lower  the  Ax^  and 
At^  values.  Fluxes  are  included  in  Eq.  D-23  along  with  the  previous 
velocity  for  linearly  approximating  the  new  velocity  V^.  In 

situations  in  which  the  dependence  is  not  close  to  linear,  it  is 
advisable  to  alter  this  dependence  accordingly.  The  consequence 
of  using  various  values  is  presented  in  Section  5 of  this  appendix. 
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0.4  MINIMIZING  COMPUTATION  TIMES 


In  that  the  computational  times  are  proportional  to  the  square 
of  the  number  of  time  intervals  used,  it  is  important  to  minimize 
the  number  of  time  intervals  used  consistent  with  accurate  predic- 
tions. In  this  section  we  shall  discuss  means  for  minimizing  the 
number  of  time  intervals. 

D . 4 . 1 Elimination  of  Early  Time  Intervals 

Two  means  were  examined  in  this  regard.  The  most  obvious 
possibility  is  to  eliminate  early  q^  when  they  cease  to  be  impor- 
tant. Unfortunately  this  is  not  particularly  fruitful  in  that 
only  rarely  do  the  effects  of  the  early  q^  become  insignificant. 

A more  practical  technique  is  to  consolidate  groups  of  the  early 
fluxes  into  single  effective  fluxes  applied  over  some  depth  during 
some  time  interval.  This  approach  has  yielded  good  results  pro- 
vided the  incident  fluxes  occur  at  the  same  depth.  However  the 
problem  of  consolidating  fluxes  applied  at  different  depths 
remains  to  be  completed  (see  Section  6 of  this  Appendix). 

D . 4 . 2 Temporary  Skipping  of  Temperature  Calculations 

Another  means  used  to  minimize  computational  times  is  based 
upon  the  fact  that  the  burning  can  become  highly  variable  over  the 
crack  elements.  In  such  situations,  much  smaller  time  steps  are 
needed  at  the  crack  locations  involving  rapid  regression  rates 
than  at  other  locations.  In  that  a single  time  Interval  is  used 
for  all  crack  elements,  there  are  instances  wherein  the  changes 
of  the  temperatures  or  regression  rates  of  some  of  the  crack 
elements  becomes  insignificant.  This  occurs  when  extremely  small 

l 

time  intervals  are  applied  to  crack  elements  involving  little  or 
no  regression. 

Whenever  such  situations  arise,  the  temperature/velocity  com- 
putations at  crack  locations  Involving  little  or  no  regression  are 
temporarily  suspended.  The  criterion  used  to  make  this  decision 
is  that  the  time  interval  required  to  produce  Important  changes 
at  the  crack  location  exceed  tha  time  Interval  being  used  by  a 
factor  of  at  least  10.  Time  intervals  temporarily  skipped  are 
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stored  for  future  use.  Temperature/velocity  computations  are 
resumed  when  the  number  of  skipped  time  intervals  exceeds  a 
specified  value  --  generally  three. 

D . 5 ACCURACY 

In  order  to  check  the  accuracy  of  the  method  we  have  applied 
it  to  a moving-boundary  problem  having  the  analytical  solution 
described  by  Eq.  8.  This  problem  involves  a constant-velocity 
boundary  that  loses  heat  convectively  to  a constant-temperature 
environment.  Initially  the  body  is  considered  to  have  a tempera- 
ture 100°C  above  its  environment.  Convective  heat  fluxes  are 
computed  using  a heat-transfer  coefficient  of  0.5  cal/cm  -sec-°C. 
Transient  fluxes  required  by  the  approximate  method  were  determined 
by  using  the  temperatures  at  the  midpoints  of  the  time  intervals 
provided  by  the  analytical  solution  described  by  Eq.  8. 

Decreases  of  temperature  obtained  by  the  two  solutions  are 
given  in  Tables  Dl,  D2  and  D3  for  velocities  of  1 , 3 and  5 cm/sec, 
respectively.  These  temperature  decreases  are  equal  to  the  dif- 
ference between  the  initial  temperature  of  100°C  minus  the  com- 
puted temperatures  within  the  body,  and  are  presented  since  they 
afford  the  most  appropriate  measure  of  the  accuracy  of  the  method. 
Time  intervals  were  constant  for  each  of  the  three  cases  treated 
and  were  determined  using  Eq.  D-23withC^  set  equal  to  0.3. 

From  the  above  tables  it  may  be  observed  that  the  agreement 
with  the  exact  temperature  changes  is  good  at  all  depths,  and 
times  for  each  of  the  velocities  considered.  Percentage  errors 
associated  with  these  results  are  shown  in  Tables  D4,  D5  and  D6. 

On  a percentage  basis  errors  are  generally  largest  at  depths  in- 
curring small  temperature  changes  of  less  than  1 C.  This  is  due 
to  taking  small  differences  between  approximations  of  the  erfc 
and  ierfc  functions.  Included,  in  the  latter  tables  in  the  con- 
sequence of  using  values  larger  than  0.3.  From  Eq.  D-23  it 
may  be  observed  that  the  value  of  0.6  corresponds  to  time  in- 
tervals four  times  larger  than  those  associated  with  a value 


D12 


SEC  WITH  C 


0 MS 

•H 


o o 

O O 

H H 

Cl  CM 

«4-  Mj 

r*v  vo 

O O 
« • 

O O 

» 4 

O O 

o a 

o o 

o o 

CM  CM 

o_o 

vo  in 
o o 

rH  O 

m 

o o 

H H 

• • 
r-»  vo 

vo  tn 

ci  'S' 

<}■  co 

o o 
• • 

oo 

rH  rH 

CM  CM 

Cl  CO 
• • 

o o 

lv  VO 

o o 

OV  00 

|v  IV 

CO  CO 

o o 

CM  CM 

Cl  Cl 

t • 

tn  m 

IV  IV 

CM  CM 

ci  ci 

Ol  00 

CJV  fv 

in  <t 

o o 

CM  CM 

tn  m 

01  Oi 

CM  Cl 

in  tn 

i H H 

H rH 

m c n 

O Ol 

VO  Sf 

Ol 

OSH 

Cl  CM 

rH  rH 

00  IV 

tntn 

CM  H 

r-  rv 

CM  H 

H H 

CM  CM 

CM  CM 

Cl  Cl 

ciO 

CM  00 

MfOO 

Cl  O 

<f  m 

Ov  Cl 

ivn. 

Cl 

00  lv 

£31  00 

c»  m 

CM  H 

CM  CM 

COCO 

Mt<t 

in  in 

VO  vO 

•MJ  VO 

HO 

CJV  Cl 

•41  C0 

tn  eo 

fv  in 

HO 

tn 

uo  in 

O O' 

O CJV 

vo  vo 

• • 

• • 

• » 

4 « 

i i 

• ) 

CO  Cl 

VO  VC 

00  00 

O CJV 

H O 

HH 

H 

H H 

H rH 

00  OQ 

VO  O 

Ov  VO 

OVVO 

VO  VO 

o oo 

MS'? 

tn<t 

H O 

fv  vO 

00  IV 

vo  tn 

i « 

I 4 

• • 

4 4 

4 4 

4 • 

oo 

in  in 

00  00 

Ov  Ol 

o a 

H H 

rH  iH 

H H 

rH  H 

HH 

CM  CM 

Cl  CM 

CM  CM 

00  VO 

<*o 

O Ol 

Cl 

01  iv 

O' 00 

m in 

CJV  ov 

*5  cm 

vtf  CM 

CM  00 

• • 

• • 

« 4 

4 • 

• 4 

1 i 

fv  OO 

cm  ci 

tn 

VO  IV 

|V  00 

00  00 

CM  CM 

Cl  CO 

C0  Cl 

COCO 

CO  CO 

Cl  Cl 

•aHsj 

iS0^ 


o m 

CM  CM 


I 


I .5 


i 

I 

I 

8 

I 


CM 

I 

O 


J3 

rt 

H 


o o 

CM  CM 

oc  00 

m m 

CO  to 

Ot  O 

■4- 

O O 

o o 

rH  rH 

CM  CH 

CM  CO 

O' 

iH 

75  00 

o o 

LT»  VO 

vo  vO 

00  00 

O'  O 

CTt  O 

a ch 

o o 

rH  rH 

CM  CM 

CO 

vj-  m 

h rv 

75  Cj  H 

a at  cm 

r H H 

co  to 

CM  CM 

O O 

r* 

CO  rH 

0)  4J  H 

“Sg  a 

!g|- 

O o.^  cr. 

o o 

fH  rH 
« » 

m co 
• » 

m m 

• t 

VO  vo 

00  00 
» » 

m in 

o o 

rH  rH 

o oo 

m cm 

to  rH 

o o 

CO  CO 

VO  VO 

OH  OO 

rH  rH 

CO  CO 

■u  ^ tj  H 

rH  rH 

rH  rH 

"»jH 

in  m 

m 

m <o 

m cm 

VO  CO 

CM  00 

qSti  °°. 

H 1 rH 

VO  VO 

rH  rH 
• • 

m m 

oo  oo 

H O 

» » 

■9  M S o 

H H 

H rH 

rH  rH 

CM  CM 

saE 

75  S <C 

m 

Mf  CM 

CM  CO 

CO  to 

o\  m 

cm  m 

rt  w 5f  vo 

xf'tf 

CO  to 

O O 

m m 

Os  O' 

CM  CM 

a ,2 

■ • 

• • 

• ■ 

• • 

« i 

• » 

M o*  oo 

rH  rH 

CM  CM 

CM  CM 

CM  CM 

to  co 

£{<*  * 
o a-g 

0 oo 

co  m 

O VO 

tO  rH 

moo 

CO  00 

Os  to 

at  H * <f 

rH  fH 

vo  m 

m m 

CM  rH 

rv  vo 

o o 

M ts 

* • 

• • 

• • 

« • 

• ♦ 

« • 

pj  1 | «3 

rH  rH 

CM  CM 

CO  CO 

m m 

rt  *J 

ra  u 

H a o 

OJ  (jS  N 

rH  VO 

rH  CO 

«fr 

OH  CM 

OO  00 

CO  iv 

o.ort  to 

uo  r> 

oo  iv 

o a> 

r~  n» 

CO  CM 

|v  to 

0 M 

• • 

• • 

• ♦ 

• • 

• » 

k « 

So  s»- 

Hg 

CM  CM 

Ht  Mf 

vo  in 

VO  vo 

IV  IV 

(v  |V 

c vo 

rv  ao 

rv  4- 

CJ>  CO 

00  00 

O'  vo 

O 00 

■H  rH 

ot  a> 

to  to 

VO  VO 

mm 

rH  H 

vo  m 

• 

• • 

• • 

< t 

k • 

t * 

• ■ 

CM 

m m 

00  00 

O'  O' 

OO 

rH  rH 

rH  rH 

rH  rH 

H H 

rH  rH 

O CM 

O'  rH 

O' 

rH  CM 

m oo 

CO  O 

OH  VO 

VO  H 

C3<t 

O'  CO 

m co 

O CO 

O 

■ • 

• 1 

t < 

t • 

• • 

1 1 

H H 

CO  t* 

m m 

m vo 

VO  vo 

iv  rv 

H rH 

rH  rH 

rH  _H_ 

rH  rH 

rH  rH 

rH  rH 

• C/3 

S | 2-0 

h 2 S 4) 

ta  a)  k2 

inon  hj  p 

in 

10 

15 

20 

25 

30 

OOfl 

H H 

mm 

O' 

00 

i-« 

VO 

m 

Mf 

ill 

o 

• 

rH 

i 

CM 

■ 

CO 

♦ 

vfr 

• 

m 

k 

wmmmmm 

!> 

»rH 

W 

u 

05 

o. 

« 

at 

u 


co 

CM 

I 

a 


CO 

rH 

I 

a 


(A 

tr 

w 


TO 

0) 

« 

<>>4 

U 

TO 

8 

rH 

o 

75 

4J 


U 

■H 


D14 


...  - .« 


00 

rH  rH 

OO 

<t  co 

r—  vo 

ON  00 

O 00 

fr 

o o 

rH  rH 

N cm 

CO  CO 

no  m 

00 

o 

CO  IN 

ON  00 

00  00 

NO  VO 

CO  O 

<o  co 

n 

CO 

O © 

r-4  H 

co  co 

m m 

n n. 

00  00 

rH 

W Oj 

A > 

4-1 

*0  ro 

O.  (U 

co 

n n» 

rH  rH 

VO  vt 

NO  CM 

CS  00 

ai  4J 
° 2 M 

els 

»gl 

M s 

00 

o o 

CO  n 

• a 

NO  vo 
• « 

00  00 
• • 

© © 
4 4 

CM  rH 

» 4 

VD 

H H 

H rH 

in 

ON  00 

no  m 

oo  in 

N 

ON  NO 

o a5* 

CO 

rH  rH 

m m 

ON  ON 

CM  CM 

n -0 

VO  VO 

S£ij 

u 8 
o dfS 

o 4)  ** 

in 

« 

m m 

oo  r* 

O vo 

i—l  i— i 

00  CO 

r-4  i— 1 
NO  r-4 

r-l  r-l 

m iH 

_ «1  M 

C 01  r 

00 

co  <o 

ON  ON 

m -o 

00  00 

r-l  rH 

coco 

« « 

• 4 

« 4 

4 • 

• * 

4 a 

•rl  U g 

r-4  r-4 

rH  rH 

CM  CM 

CM  CS 

a $ 

S«2* 

W 2 

o 

CO  r— 1 

no 

Ntf  O 

o m 

rH  r*H 

<f  O 

« « sr  o> 

n.  p> 

NO  NO 

CM  CS 

r-  vo 

O ON 

CM  CM 

E S* 
“§•1 

CO 

■ » 

H H 

4 4 

CM  CS 

CM  N 

co  CM 

COCO 

CO 

r-4  ON 

r-  r- 

OO  CO 

o o 

NO  VO 

CM  r-4 

04  H ^ 

ON 

-t  co 

in  m 

CM  CM 

00  00 

r-4  r-* 

<r  vt 

M a 

* 

• • 

• • 

4 • 

4 • 

4 4 

4 a 

3 ' E 

<N 

r-4  rH 

CM  CM 

coco 

CO  CO 

•d-Mf 

-os- 

rt  ' s 
^ « c 

in 

oo  m 

CM  VO 

vo  m 

CM  00 

H P** 

om 

on 

m m 

© ON 

00  00 

Mf  CO 

oo  r-> 

rH  O 

« 

• • 

• • 

4 • 

• 4 

4 4 

• 4 

a o 

rH 

N CM 

<r  co 

Sf-  sc 

min 

n m 

vo  >o 

d 

00 

rIN 

<s  m 

co  n- 

m co 

in 

VO  <f 

•H 

ON 

MfMt 

o © 

On  On 

m in 

ON  ON 

CM  CM 

A 

• • 

i i 

1 • 

■ » 

4 4 

• 4 

sfsf 

NO  VO 

vo  VO 

P--  N 

p*.  n« 

00  00 

o o> 

NO  N 

N K 

ON  <t 

CO  CO 

vo  <t 

O N 

NO  ON 

in  oo 

r-4  3 

NO  00 

ON  H 

o 

4 4 

• 1 

• • 

4 4 

4 4 

i*n  n. 

CO  00 

On  ON 

© d 

OO 

© H 

"™n  r 

i-l  H 

H rH 

rH  H 

T^™5 

mmmam 

mmtmm 

mmmmm 

0 OJr- 

Mg 

'd 

•H  P 

0) 

n 

o 

m 

o 

m 

© 

£ 

rH 

H 

CM 

CM 

CO 

OOP 

H h- 

•0 

01 

A 

w d 

0 

co 

m 

oo 

O 

in 

m 

a l 

04 

w 

n 

vo 

ON 

3 

NO 

ON 

rt  *h 

o 

o 

o 

H 

r-4 

tH  h 

a 

■ 

1 

• 

i 

a 

4 

W 

-t-  * 


f0 

P.  o o 

o 

H- 

O O 

0 

CO 

P P 

Q 

*d 

CO  CO 

a 

h*  rt  rt 

a 

Pi 

6>  p 

n> 

O 

p p 

p 

£5 

rt  rt 

P. 

9 

(t> 

0 

rt  rt 

P* 

P 

P*  P* 

rt 

Pi  Pi 

< 

rt  rt 

P 

> 

H* 

X 

(0  (0 

c 

H‘ 

CO  CO 

n> 

« 

rt  rt 

» 6»  pi 

0.0"  cr 

Ml  0 M4  H 
O tt  H'  H' 
H ft  CO  CO 

erp*  cr* 

O co  <0  (D 
co  p. 


> 

o x p.  « 

M M>  (0  H> 
►O  N 
rt  (D 

tr 

co  a 


rt 

Z E’ 

rt 

hi  ro 

ff  3 


(0 

H 

CO 

CO 

w 

»p| 

H* 

CO 

CD 

61 

(0 

9 

n 

Eq 

61 

i 

*P 

*P 

o 

P 

i 

N» 

CO 

CO 

o. 

*0  CVOJ  4>  ■t*  CO  ts>  I— 1 
S3  4>  N3  00  O to  -p.  Ov  00 


M UNSMHH 
Cn  OCn  OUlOWOUl 


<J\  U>  U>  CO  U)  OJ  NJ  U1  IsJ 
Oi  00  'O  O^H^OOUl 


UJ  CO  4>  CO  4>  CO 
• ••••• 
VO  ^ >0  to  kO 


UlOUTlffiNJ-O 

4>>VflVOClHM 


I I I 
I l | hhn 
• *•••» 

oo  oo  vo  cn  cn  m 


Ul  I I I 


Cn  I cncol  MHHfi'i'Vfl 


to  O'  00  hnmhnn 

• I • «»••*« 

CO  00  O VO  O O 00  Q CO 


w 

0 HH 

co  h-  p 

co  g "p 

O (D  w 
» (t> 

o. 


H 1-3 
9 0 0 
rt  Hi  rt 

» 2 _>  p. 
i*32H 
o-S  g sa 
t-*n>  o 

C0 


N>  0 0$ 
9 o 
O M 
<3  9 9 
61  CD  9 
9 PJ  O 

H«  to  H 

0 ffi  <0 
O C « 

• CO  O 
CO  H*  Ml 

<9*6 
6i  o re 
HO  9 
O d p 6> 

• S rt  rt 
cn  co  d 


6l  v 9 
9 § to 
H»  2. 

U! 


s-S  " 

gss 

CO 


F>‘. 


■ '«**#« 


Table  D-4 


Table  D-5 


oo  intsi  n p-(  cs 
CS  CS  CS  CS  CS  CS  00 


OOH  rt  »4  r-OH 
till 


in^nNOW 

NHHrIH  I 
I t I I 1 


ts  o>  cl  cs  n.  vo 
oo  in  co  on  cs  cs 


oornHsi-fv  o 

ClrtrlrtHH  CO  CM  ID 


oo  m v© co»n oo  ho  m 

**•»•«  ' ■ I • 

NNHHiHiH  COOS 


m © m o in  o 

rH  rH  CS  CS  CO 


0>00h*v©m^  CO  CS  r-t 

©p-uscn^m  cnrs  oo 


T3 

0) 

on  *rl 
CS  H 
• Pu 
Q fM 
(0 


cr 

0) 

w 

M 

cd 

« 

o> 

•H 

in 

1 cr 

* 

w 

m 

flJ 

rH 

g 

9 

rS 

M 

4-1 

01 

W 

01 

B 

X 

•H 

V 

Ql 

X 

jj 

0 -rl 
■rj  j*j 

u 

^ < 

> CS 

0) 

o 

X) 

u 

u 

nt  II 

4-1 

o)  m 

O 

XXI 

V AJ 

ou 

Or  Or 

H 

0)  01  CS 

•r( 

uu  o 

M , 

HI  £ «rl  C 
£ UH  «l 
4J  O 

M <tf  4J  o" 
« jdg  , 

Xt  a o & M 

M -rl  E O 
■rl  i-H  <U  44 
rl  J O 
,o  ofl  rt  w 

0 «rl  0) 

w Mft  3 

W # M rl 
0)  -r4  <0 

| ^ 

Hi! 

§ s3  8 

CJ  U S*  (2 

* 4-  ? 


* 


I 


r 


in 

o in  i-4  sj  r--  ms 

i 

tfi 

• 

tl  44  «}  ° 

f4  D 

& Sd 

CO  CM  CM  CM  CM  CM 

I 

CM 

rH 

o ^ u n 

mt  cm  co  cm  n n 

^3 

V0 

|J  rj  * 

^ .5  w o 

1 1 1 

1 1* 

CO 

(3  H 3 

1 

•rl  u O 
e 'H 

”S» 

8 g ^ 

E Vl  CM 

W u 0 

cm  n vf  os  vo  m 

O 00 

o 

44  O 

CM  rH  r-4  1 1 1 

in  co 

r-H 

1 1 1 

i i 

1 

04 

V 

a w in 

OS  Sf  VO  O 00 

o 

44  44  OJ 

• ••••• 

* 

nj  a 3 o 

o-  «o-  CO  CO  CO 

00 

M JIH 
0)  o to 

Oh  44  > 

Sfio-  4 

r-4 

r-J 

H CM* 

fio  * 

44  »h  co 

oo  scMinoa 

in  so 

OS 

O M 

« • 

• 

(0  9 0 
(0  0)  o 

M (0  -H 
O (0.  Vi 
Vt  0)  <8 

CO  r-4  iH  i-4  CM  r-J 

CO  CM 

Vi  V4  > 

W o 

0J  Vi 

COO  CM 

00  CSCOOMf  CM 

CM  0- 

CO 

n 4-i 

• t « • • • 

• * 

f 

£ ° 

N H H H i~<  H 

co  cM 

00 

• n 

■e 

0 OJr-4 

S-o 

i 

«H  S£ 

44  4-1  44  B 

n o in  o n o 

n 

H H CM  CM  CO 

n 

0 0 c 

H H 

■ 

T3 

OJ  * 
to  OJ  O 

co  uo  o co  in 

o © 

co 

o.  a oj 

co  so  os  co  so  Oi 

CO  VO 

Os 

|J'3« 

OOOHHH 

r-4  CM 

CM 

SJH§ 

• • 

* , 

1 

VO 

OS 

H 

» 

• 

• 

(J 

o 

o 

o 

D18 


0 MM 
43  43 

01  44  4J 

N ex  a. 

•H  0>  4)  CM 

W tJ  tJ  U 

o)  o)  » 

43  43  0 

44  44  M tfl 

< 

W CO  i-4 

41  01  44  U 

43  43  ?3 

(0  (0  0)  M 

•r4  t4  a o 
r4  r— I 3 *44 

"rt  (fl  M 

44  4->  H 4> 

w w a*  a 

04  * 

44  44  -0  > 


I 

! 


I 


i 

I 


I 


r wi^qrg; 


of  0.3  while  the  value  of  Q.9  increases  the  time  intervals  by 
a factor  of  nine.  As  one  would  expect,  increased  values  usually 
increase  the  errors.  Also  it  should  be  observed  that  the  effect 
of  varying  C2  is  most  pronounced  with  the  larger  values.  Thi3 
is  not  at  all  surprising  in  that  lower  values  yield  smaller  dis- 

placements Ax^.  The  result  is  a lessening  of  the  thermal  resistance 
across  the  displacements,  and  hence  of  the  importance  of  where  the 
source  fluxes  q^  are  located  within  Ax^. 

In  view  of  the  above  results,  it  is  advisable  to  use  values 
of  0.3  or  less  depending  upon  the  accuracy  desired.  Here  one  should 
be  particularly  cognizant  of  trade-offs  between  accuracy  and  com- 
puter execution  times  in  that  the  execution  times  associated  with 
a given  period  of  time  are  roughly  proportional  to  1/C^.  Means 
for  speeding  the  execution  times  are  discussed  in  the  next  section. 


D. 6 SUMMARY 

In  summary,  the  time- dependent  flux  is  approximated  by  a 
series  of  constant  fluxes  applied  during  various  time  inter- 
vals At^,  In  order  to  use  solutions  for  a stationary  semi-infinite 
body,  it  is  necessary  to  account  for  two  facts.  The  first  is  that 
there  should  be  no  fluxes  conducted  across  the  displaced  boundary 
from  previous  fluxes.  The  second  is  that  in  .depth  heating  caused 
by  prior  fluxes  must  be  maintained.  To  achieve  these  objectives, 
an  additional  flux  is  introduced  to  cancel  the  undesirable  heat 
conducted,  across  a given  depth  from  earlier  fluxes  during  the  most 
recent  time  interval.  This  depth  corresponds  to  where  the  next 
constant  flux  is  to  be  applied  and  will  be  discussed  in  the  next 
paragraph.  The  result  is  a flux  lower  than  q^  by  the  mean  con- 
ductive flux  arising  from  prior  q^. 

Next  it  is  necessary  to  approximate  the  heat  conducted  across 
the  depth  at  which  the  boundary  ultimately  arrives  at  the  end  of 
the  time  interval.  This  is  accomplished  by  selecting  sufficiently 
small  time  intervals  of  values  so  that  the  thermal  resistance 
across  displacement  Ax^  is  small.  By  doing  so  the  variations  of 
the  conductive  fluxes  across  the  displacements  are  also  kept 
small.  Fluxes  q^  are  then  located  at  an  effective  depth  within 
the  small  displacements  according  to  the  value  assigned  the 
constant  C2 • 

Finally  Duhamel's  principle  is  UBed  to  account  for  the  fact 
that  each  of  the  fluxes  exist  only  over  given  time  intervals 
At^.  By  this  technique  one  arrives  at  a series  of  constant  fluxes 
q^  applied  at  various  depths  x^  and  time  intervals  Atj.  with  which 
to  determine  the  temperature  changes  of  the  regressing  propellant. 

The  advantage  of  thiB  method  over  finite-difference  techniques 
is  threefold.  The  first  is  that  temperature  predictions  require 
less  stored  information  than  needed  by  finite-difference  methods. 

This  is  Important  In  treating  numerous  crack  elements.  Secondly, 
the  method  provides  reasonably  accurate  predictions  of  the  temperatures 

D20 


using  relatively  large  time  intervals.  Finally,  the  method  does 
not  require  periodic  refinements  of  the  spacial  grid  to  accomo- 
date pronounced  changes  in  the  thermal  gradients  with  depth  as 
the  burning  accelerates. 

The  most  undesirable  feature  of  the  method  is  that  the  compu- 
ter execution  times  increase  approximately  linearly  with  the  num- 
ber of  time  intervals  involved  in  each  set  of  computations.  Means 
have  been  formulated  for  limiting  the  number  of  time  intervals  by 
first  determining  the  amount  of  heat  conducted  across  the  most 
recent  depth  x^  during  each  time  interval.  Then  the  time -dependent 
heat  conduction  is  approximated  by  a series  of  constant  fluxes  over 
one  or  more  time  intervals.  During  early  times,  wherein  the  fluxes 
are  least  important,  the  time  intervals  are  consolidated  yielding 
a courser  approximation  of  the  fljux  at  the  early  times.  Prelim- 
-iD.ajJr'l  computer  runs  suggest  it  is  possible  to  limit  the  number 
of  time  intervals  being  used  to  about  30  without  introducing  errors 
of  more  than  one  percent  in  the  temperature  rises.  However,  cri- 
teria for  approximating  the  fluutes  at  early  times  in  terms  of 
various  dynamic  burn  conditions  has  yet  to  be  established. 
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This  study  Involves  the  development  of  s coe^uterlsed  model  with  which  to 
identify  the  cautes  of  solid  propellant  rocket  motor  detonations  during  firing. 
The  model  treats  the  problert  of  suddenly  exposing  propellant  cracks  (or  flaws  or 
debonda)  to  chambor  gases  at  higher  pressures  and  temperatures  than  gasea 
Initially  In  ths  crack.  Simplified  analyses  are  used  to  compute  the  transient 
pressures  produced  by  the  gas  flows,  burning  and  crack  deformations  (positive 
and  negative).  These  prodictions  are  made  in  terms  of  a variety  of  parameters, 
associated  with  the  crack,  chamber,  motor  case  and  propellant.  Two  phenomena 
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appear  important  in  producing  pronounced  pressure  rises.  These  are:  (1)  partial 

complete  crack  collapse  due  to  stress  waves  reflected*  from  "the  motor  case  and 
(2)  rapid  acceleration  of  the  burning  during  the  final  stages  of  crack  collapse, 
fhe  latter  is  due  to  the  relatively  large  quantities  of  heat  acquired  by  the 
•ropel  I ant  during  the  relatively  long  periods  of  gas  filling.  Quantification  of 
*h#eu  effects  remains  to  be  accomplished. 
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